home *** CD-ROM | disk | FTP | other *** search
/ EnigmA Amiga Run 1995 November / EnigmA AMIGA RUN 02 (1995)(G.R. Edizioni)(IT)[!][issue 1995-11][Skylink CD].iso / earcd / program / gcc / gcc270-s.lha / gcc-2.7.0-amiga / real.c < prev    next >
C/C++ Source or Header  |  1995-08-24  |  128KB  |  6,121 lines
  1. /* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
  2.    and support for XFmode IEEE extended real floating point arithmetic.
  3.    Copyright (C) 1993, 1994, 1995 Free Software Foundation, Inc.
  4.    Contributed by Stephen L. Moshier (moshier@world.std.com).
  5.  
  6. This file is part of GNU CC.
  7.  
  8. GNU CC is free software; you can redistribute it and/or modify
  9. it under the terms of the GNU General Public License as published by
  10. the Free Software Foundation; either version 2, or (at your option)
  11. any later version.
  12.  
  13. GNU CC is distributed in the hope that it will be useful,
  14. but WITHOUT ANY WARRANTY; without even the implied warranty of
  15. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16. GNU General Public License for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with GNU CC; see the file COPYING.  If not, write to
  20. the Free Software Foundation, 59 Temple Place - Suite 330,
  21. Boston, MA 02111-1307, USA.  */
  22.  
  23. #include <stdio.h>
  24. #include <errno.h>
  25. #include "config.h"
  26. #include "tree.h"
  27.  
  28. #ifndef errno
  29. extern int errno;
  30. #endif
  31.  
  32. /* To enable support of XFmode extended real floating point, define
  33. LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
  34.  
  35. To support cross compilation between IEEE, VAX and IBM floating
  36. point formats, define REAL_ARITHMETIC in the tm.h file.
  37.  
  38. In either case the machine files (tm.h) must not contain any code
  39. that tries to use host floating point arithmetic to convert
  40. REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
  41. etc.  In cross-compile situations a REAL_VALUE_TYPE may not
  42. be intelligible to the host computer's native arithmetic.
  43.  
  44. The emulator defaults to the host's floating point format so that
  45. its decimal conversion functions can be used if desired (see
  46. real.h).
  47.  
  48. The first part of this file interfaces gcc to a floating point
  49. arithmetic suite that was not written with gcc in mind.  Avoid
  50. changing the low-level arithmetic routines unless you have suitable
  51. test programs available.  A special version of the PARANOIA floating
  52. point arithmetic tester, modified for this purpose, can be found on
  53. usc.edu: /pub/C-numanal/ieeetest.zoo.  Other tests, and libraries of
  54. XFmode and TFmode transcendental functions, can be obtained by ftp from
  55. netlib.att.com: netlib/cephes.   */
  56.  
  57. /* Type of computer arithmetic.
  58.    Only one of DEC, IBM, IEEE, or UNK should get defined.
  59.  
  60.    `IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
  61.    to big-endian IEEE floating-point data structure.  This definition
  62.    should work in SFmode `float' type and DFmode `double' type on
  63.    virtually all big-endian IEEE machines.  If LONG_DOUBLE_TYPE_SIZE
  64.    has been defined to be 96, then IEEE also invokes the particular
  65.    XFmode (`long double' type) data structure used by the Motorola
  66.    680x0 series processors.
  67.  
  68.    `IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
  69.    little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
  70.    has been defined to be 96, then IEEE also invokes the particular
  71.    XFmode `long double' data structure used by the Intel 80x86 series
  72.    processors.
  73.  
  74.    `DEC' refers specifically to the Digital Equipment Corp PDP-11
  75.    and VAX floating point data structure.  This model currently
  76.    supports no type wider than DFmode.
  77.  
  78.    `IBM' refers specifically to the IBM System/370 and compatible
  79.    floating point data structure.  This model currently supports
  80.    no type wider than DFmode.  The IBM conversions were contributed by
  81.    frank@atom.ansto.gov.au (Frank Crawford).
  82.  
  83.    If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
  84.    then `long double' and `double' are both implemented, but they
  85.    both mean DFmode.  In this case, the software floating-point
  86.    support available here is activated by writing
  87.       #define REAL_ARITHMETIC
  88.    in tm.h. 
  89.  
  90.    The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
  91.    and may deactivate XFmode since `long double' is used to refer
  92.    to both modes.
  93.  
  94.    The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
  95.    contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
  96.    separate the floating point unit's endian-ness from that of
  97.    the integer addressing.  This permits one to define a big-endian
  98.    FPU on a little-endian machine (e.g., ARM).  An extension to
  99.    BYTES_BIG_ENDIAN may be required for some machines in the future.
  100.    These optional macros may be defined in tm.h.  In real.h, they
  101.    default to WORDS_BIG_ENDIAN, etc., so there is no need to define
  102.    them for any normal host or target machine on which the floats
  103.    and the integers have the same endian-ness.   */
  104.  
  105.  
  106. /* The following converts gcc macros into the ones used by this file.  */
  107.  
  108. /* REAL_ARITHMETIC defined means that macros in real.h are
  109.    defined to call emulator functions.  */
  110. #ifdef REAL_ARITHMETIC
  111.  
  112. #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  113. /* PDP-11, Pro350, VAX: */
  114. #define DEC 1
  115. #else /* it's not VAX */
  116. #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
  117. /* IBM System/370 style */
  118. #define IBM 1
  119. #else /* it's also not an IBM */
  120. #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
  121. #define IEEE
  122. #else /* it's not IEEE either */
  123. /* UNKnown arithmetic.  We don't support this and can't go on. */
  124. unknown arithmetic type
  125. #define UNK 1
  126. #endif /* not IEEE */
  127. #endif /* not IBM */
  128. #endif /* not VAX */
  129.  
  130. #define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
  131.  
  132. #else
  133. /* REAL_ARITHMETIC not defined means that the *host's* data
  134.    structure will be used.  It may differ by endian-ness from the
  135.    target machine's structure and will get its ends swapped
  136.    accordingly (but not here).  Probably only the decimal <-> binary
  137.    functions in this file will actually be used in this case.  */
  138.  
  139. #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  140. #define DEC 1
  141. #else /* it's not VAX */
  142. #if HOST_FLOAT_FORMAT == IBM_FLOAT_FORMAT
  143. /* IBM System/370 style */
  144. #define IBM 1
  145. #else /* it's also not an IBM */
  146. #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
  147. #define IEEE
  148. #else /* it's not IEEE either */
  149. unknown arithmetic type
  150. #define UNK 1
  151. #endif /* not IEEE */
  152. #endif /* not IBM */
  153. #endif /* not VAX */
  154.  
  155. #define REAL_WORDS_BIG_ENDIAN HOST_FLOAT_WORDS_BIG_ENDIAN
  156.  
  157. #endif /* REAL_ARITHMETIC not defined */
  158.  
  159. /* Define INFINITY for support of infinity.
  160.    Define NANS for support of Not-a-Number's (NaN's).  */
  161. #if !defined(DEC) && !defined(IBM)
  162. #define INFINITY
  163. #define NANS
  164. #endif
  165.  
  166. /* Support of NaNs requires support of infinity. */
  167. #ifdef NANS
  168. #ifndef INFINITY
  169. #define INFINITY
  170. #endif
  171. #endif
  172.  
  173. /* Find a host integer type that is at least 16 bits wide,
  174.    and another type at least twice whatever that size is. */
  175.  
  176. #if HOST_BITS_PER_CHAR >= 16
  177. #define EMUSHORT char
  178. #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
  179. #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
  180. #else
  181. #if HOST_BITS_PER_SHORT >= 16
  182. #define EMUSHORT short
  183. #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
  184. #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
  185. #else
  186. #if HOST_BITS_PER_INT >= 16
  187. #define EMUSHORT int
  188. #define EMUSHORT_SIZE HOST_BITS_PER_INT
  189. #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
  190. #else
  191. #if HOST_BITS_PER_LONG >= 16
  192. #define EMUSHORT long
  193. #define EMUSHORT_SIZE HOST_BITS_PER_LONG
  194. #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
  195. #else
  196. /*  You will have to modify this program to have a smaller unit size. */
  197. #define EMU_NON_COMPILE
  198. #endif
  199. #endif
  200. #endif
  201. #endif
  202.  
  203. #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
  204. #define EMULONG short
  205. #else
  206. #if HOST_BITS_PER_INT >= EMULONG_SIZE
  207. #define EMULONG int
  208. #else
  209. #if HOST_BITS_PER_LONG >= EMULONG_SIZE
  210. #define EMULONG long
  211. #else
  212. #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
  213. #define EMULONG long long int
  214. #else
  215. /*  You will have to modify this program to have a smaller unit size. */
  216. #define EMU_NON_COMPILE
  217. #endif
  218. #endif
  219. #endif
  220. #endif
  221.  
  222.  
  223. /* The host interface doesn't work if no 16-bit size exists. */
  224. #if EMUSHORT_SIZE != 16
  225. #define EMU_NON_COMPILE
  226. #endif
  227.  
  228. /* OK to continue compilation. */
  229. #ifndef EMU_NON_COMPILE
  230.  
  231. /* Construct macros to translate between REAL_VALUE_TYPE and e type.
  232.    In GET_REAL and PUT_REAL, r and e are pointers.
  233.    A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
  234.    in memory, with no holes.  */
  235.  
  236. #if LONG_DOUBLE_TYPE_SIZE == 96
  237. /* Number of 16 bit words in external e type format */
  238. #define NE 6
  239. #define MAXDECEXP 4932
  240. #define MINDECEXP -4956
  241. #define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
  242. #define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
  243. #else /* no XFmode */
  244. #if LONG_DOUBLE_TYPE_SIZE == 128
  245. #define NE 10
  246. #define MAXDECEXP 4932
  247. #define MINDECEXP -4977
  248. #define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
  249. #define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
  250. #else
  251. #define NE 6
  252. #define MAXDECEXP 4932
  253. #define MINDECEXP -4956
  254. #ifdef REAL_ARITHMETIC
  255. /* Emulator uses target format internally
  256.    but host stores it in host endian-ness. */
  257.  
  258. #define GET_REAL(r,e)                        \
  259. do {                                \
  260.      if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN)    \
  261.        e53toe ((unsigned EMUSHORT*) (r), (e));            \
  262.      else                            \
  263.        {                            \
  264.      unsigned EMUSHORT w[4];                \
  265.      w[3] = ((EMUSHORT *) r)[0];                \
  266.      w[2] = ((EMUSHORT *) r)[1];                \
  267.      w[1] = ((EMUSHORT *) r)[2];                \
  268.      w[0] = ((EMUSHORT *) r)[3];                \
  269.      e53toe (w, (e));                    \
  270.        }                            \
  271.    } while (0)
  272.  
  273. #define PUT_REAL(e,r)                        \
  274. do {                                \
  275.      if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN)    \
  276.        etoe53 ((e), (unsigned EMUSHORT *) (r));            \
  277.      else                            \
  278.        {                            \
  279.      unsigned EMUSHORT w[4];                \
  280.      etoe53 ((e), w);                    \
  281.      *((EMUSHORT *) r) = w[3];                \
  282.      *((EMUSHORT *) r + 1) = w[2];                \
  283.      *((EMUSHORT *) r + 2) = w[1];                \
  284.      *((EMUSHORT *) r + 3) = w[0];                \
  285.        }                            \
  286.    } while (0)
  287.  
  288. #else /* not REAL_ARITHMETIC */
  289.  
  290. /* emulator uses host format */
  291. #define GET_REAL(r,e) e53toe ((unsigned EMUSHORT *) (r), (e))
  292. #define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))
  293.  
  294. #endif /* not REAL_ARITHMETIC */
  295. #endif /* not TFmode */
  296. #endif /* no XFmode */
  297.  
  298.  
  299. /* Number of 16 bit words in internal format */
  300. #define NI (NE+3)
  301.  
  302. /* Array offset to exponent */
  303. #define E 1
  304.  
  305. /* Array offset to high guard word */
  306. #define M 2
  307.  
  308. /* Number of bits of precision */
  309. #define NBITS ((NI-4)*16)
  310.  
  311. /* Maximum number of decimal digits in ASCII conversion
  312.  * = NBITS*log10(2)
  313.  */
  314. #define NDEC (NBITS*8/27)
  315.  
  316. /* The exponent of 1.0 */
  317. #define EXONE (0x3fff)
  318.  
  319. extern int extra_warnings;
  320. extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
  321. extern unsigned EMUSHORT elog2[], esqrt2[];
  322.  
  323. static void endian    PROTO((unsigned EMUSHORT *, long *,
  324.                    enum machine_mode));
  325. static void eclear    PROTO((unsigned EMUSHORT *));
  326. static void emov    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  327. static void eabs    PROTO((unsigned EMUSHORT *));
  328. static void eneg    PROTO((unsigned EMUSHORT *));
  329. static int eisneg    PROTO((unsigned EMUSHORT *));
  330. static int eisinf    PROTO((unsigned EMUSHORT *));
  331. static int eisnan    PROTO((unsigned EMUSHORT *));
  332. static void einfin    PROTO((unsigned EMUSHORT *));
  333. static void enan    PROTO((unsigned EMUSHORT *, int));
  334. static void emovi    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  335. static void emovo    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  336. static void ecleaz    PROTO((unsigned EMUSHORT *));
  337. static void ecleazs    PROTO((unsigned EMUSHORT *));
  338. static void emovz    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  339. static void einan    PROTO((unsigned EMUSHORT *));
  340. static int eiisnan    PROTO((unsigned EMUSHORT *));
  341. static int eiisneg    PROTO((unsigned EMUSHORT *));
  342. static void eiinfin    PROTO((unsigned EMUSHORT *));
  343. static int eiisinf    PROTO((unsigned EMUSHORT *));
  344. static int ecmpm    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  345. static void eshdn1    PROTO((unsigned EMUSHORT *));
  346. static void eshup1    PROTO((unsigned EMUSHORT *));
  347. static void eshdn8    PROTO((unsigned EMUSHORT *));
  348. static void eshup8    PROTO((unsigned EMUSHORT *));
  349. static void eshup6    PROTO((unsigned EMUSHORT *));
  350. static void eshdn6    PROTO((unsigned EMUSHORT *));
  351. static void eaddm    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  352. static void esubm    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  353. static void m16m    PROTO((unsigned int, unsigned short *,
  354.                    unsigned short *));
  355. static int edivm    PROTO((unsigned short *, unsigned short *));
  356. static int emulm    PROTO((unsigned short *, unsigned short *));
  357. static void emdnorm    PROTO((unsigned EMUSHORT *, int, int, EMULONG, int));
  358. static void esub    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  359.                    unsigned EMUSHORT *));
  360. static void eadd    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  361.                    unsigned EMUSHORT *));
  362. static void eadd1    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  363.                    unsigned EMUSHORT *));
  364. static void ediv    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  365.                    unsigned EMUSHORT *));
  366. static void emul    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  367.                    unsigned EMUSHORT *));
  368. static void e53toe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  369. static void e64toe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  370. static void e113toe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  371. static void e24toe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  372. static void etoe113    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  373. static void toe113    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  374. static void etoe64    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  375. static void toe64    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  376. static void etoe53    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  377. static void toe53    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  378. static void etoe24    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  379. static void toe24    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  380. static int ecmp        PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  381. static void eround    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  382. static void ltoe    PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *));
  383. static void ultoe    PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *));
  384. static void eifrac    PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *,
  385.                    unsigned EMUSHORT *));
  386. static void euifrac    PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *,
  387.                    unsigned EMUSHORT *));
  388. static int eshift    PROTO((unsigned EMUSHORT *, int));
  389. static int enormlz    PROTO((unsigned EMUSHORT *));
  390. static void e24toasc    PROTO((unsigned EMUSHORT *, char *, int));
  391. static void e53toasc    PROTO((unsigned EMUSHORT *, char *, int));
  392. static void e64toasc    PROTO((unsigned EMUSHORT *, char *, int));
  393. static void e113toasc    PROTO((unsigned EMUSHORT *, char *, int));
  394. static void etoasc    PROTO((unsigned EMUSHORT *, char *, int));
  395. static void asctoe24    PROTO((char *, unsigned EMUSHORT *));
  396. static void asctoe53    PROTO((char *, unsigned EMUSHORT *));
  397. static void asctoe64    PROTO((char *, unsigned EMUSHORT *));
  398. static void asctoe113    PROTO((char *, unsigned EMUSHORT *));
  399. static void asctoe    PROTO((char *, unsigned EMUSHORT *));
  400. static void asctoeg    PROTO((char *, unsigned EMUSHORT *, int));
  401. static void efloor    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  402. static void efrexp    PROTO((unsigned EMUSHORT *, int *,
  403.                    unsigned EMUSHORT *));
  404. static void eldexp    PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *));
  405. static void eremain    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  406.                    unsigned EMUSHORT *));
  407. static void eiremain    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  408. static void mtherr    PROTO((char *, int));
  409. static void dectoe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  410. static void etodec    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  411. static void todec    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  412. static void ibmtoe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  413.                    enum machine_mode));
  414. static void etoibm    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  415.                    enum machine_mode));
  416. static void toibm    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  417.                    enum machine_mode));
  418. static void make_nan    PROTO((unsigned EMUSHORT *, int, enum machine_mode));
  419. static void uditoe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  420. static void ditoe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  421. static void etoudi    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  422. static void etodi    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  423. static void esqrt    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  424.  
  425. /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
  426.    swapping ends if required, into output array of longs.  The
  427.    result is normally passed to fprintf by the ASM_OUTPUT_ macros.   */
  428.  
  429. static void 
  430. endian (e, x, mode)
  431.      unsigned EMUSHORT e[];
  432.      long x[];
  433.      enum machine_mode mode;
  434. {
  435.   unsigned long th, t;
  436.  
  437.   if (REAL_WORDS_BIG_ENDIAN)
  438.     {
  439.       switch (mode)
  440.     {
  441.  
  442.     case TFmode:
  443.       /* Swap halfwords in the fourth long. */
  444.       th = (unsigned long) e[6] & 0xffff;
  445.       t = (unsigned long) e[7] & 0xffff;
  446.       t |= th << 16;
  447.       x[3] = (long) t;
  448.  
  449.     case XFmode:
  450.  
  451.       /* Swap halfwords in the third long. */
  452.       th = (unsigned long) e[4] & 0xffff;
  453.       t = (unsigned long) e[5] & 0xffff;
  454.       t |= th << 16;
  455.       x[2] = (long) t;
  456.       /* fall into the double case */
  457.  
  458.     case DFmode:
  459.  
  460.       /* swap halfwords in the second word */
  461.       th = (unsigned long) e[2] & 0xffff;
  462.       t = (unsigned long) e[3] & 0xffff;
  463.       t |= th << 16;
  464.       x[1] = (long) t;
  465.       /* fall into the float case */
  466.  
  467.     case HFmode:
  468.     case SFmode:
  469.  
  470.       /* swap halfwords in the first word */
  471.       th = (unsigned long) e[0] & 0xffff;
  472.       t = (unsigned long) e[1] & 0xffff;
  473.       t |= th << 16;
  474.       x[0] = t;
  475.       break;
  476.  
  477.     default:
  478.       abort ();
  479.     }
  480.     }
  481.   else
  482.     {
  483.       /* Pack the output array without swapping. */
  484.  
  485.       switch (mode)
  486.     {
  487.  
  488.     case TFmode:
  489.  
  490.       /* Pack the fourth long. */
  491.       th = (unsigned long) e[7] & 0xffff;
  492.       t = (unsigned long) e[6] & 0xffff;
  493.       t |= th << 16;
  494.       x[3] = (long) t;
  495.  
  496.     case XFmode:
  497.  
  498.       /* Pack the third long.
  499.          Each element of the input REAL_VALUE_TYPE array has 16 useful bits
  500.          in it.  */
  501.       th = (unsigned long) e[5] & 0xffff;
  502.       t = (unsigned long) e[4] & 0xffff;
  503.       t |= th << 16;
  504.       x[2] = (long) t;
  505.       /* fall into the double case */
  506.  
  507.     case DFmode:
  508.  
  509.       /* pack the second long */
  510.       th = (unsigned long) e[3] & 0xffff;
  511.       t = (unsigned long) e[2] & 0xffff;
  512.       t |= th << 16;
  513.       x[1] = (long) t;
  514.       /* fall into the float case */
  515.  
  516.     case HFmode:
  517.     case SFmode:
  518.  
  519.       /* pack the first long */
  520.       th = (unsigned long) e[1] & 0xffff;
  521.       t = (unsigned long) e[0] & 0xffff;
  522.       t |= th << 16;
  523.       x[0] = t;
  524.       break;
  525.  
  526.     default:
  527.       abort ();
  528.     }
  529.     }
  530. }
  531.  
  532.  
  533. /* This is the implementation of the REAL_ARITHMETIC macro.  */
  534.  
  535. void 
  536. earith (value, icode, r1, r2)
  537.      REAL_VALUE_TYPE *value;
  538.      int icode;
  539.      REAL_VALUE_TYPE *r1;
  540.      REAL_VALUE_TYPE *r2;
  541. {
  542.   unsigned EMUSHORT d1[NE], d2[NE], v[NE];
  543.   enum tree_code code;
  544.  
  545.   GET_REAL (r1, d1);
  546.   GET_REAL (r2, d2);
  547. #ifdef NANS
  548. /*  Return NaN input back to the caller. */
  549.   if (eisnan (d1))
  550.     {
  551.       PUT_REAL (d1, value);
  552.       return;
  553.     }
  554.   if (eisnan (d2))
  555.     {
  556.       PUT_REAL (d2, value);
  557.       return;
  558.     }
  559. #endif
  560.   code = (enum tree_code) icode;
  561.   switch (code)
  562.     {
  563.     case PLUS_EXPR:
  564.       eadd (d2, d1, v);
  565.       break;
  566.  
  567.     case MINUS_EXPR:
  568.       esub (d2, d1, v);        /* d1 - d2 */
  569.       break;
  570.  
  571.     case MULT_EXPR:
  572.       emul (d2, d1, v);
  573.       break;
  574.  
  575.     case RDIV_EXPR:
  576. #ifndef REAL_INFINITY
  577.       if (ecmp (d2, ezero) == 0)
  578.     {
  579. #ifdef NANS
  580.     enan (v, eisneg (d1) ^ eisneg (d2));
  581.     break;
  582. #else
  583.     abort ();
  584. #endif
  585.     }
  586. #endif
  587.       ediv (d2, d1, v);    /* d1/d2 */
  588.       break;
  589.  
  590.     case MIN_EXPR:        /* min (d1,d2) */
  591.       if (ecmp (d1, d2) < 0)
  592.     emov (d1, v);
  593.       else
  594.     emov (d2, v);
  595.       break;
  596.  
  597.     case MAX_EXPR:        /* max (d1,d2) */
  598.       if (ecmp (d1, d2) > 0)
  599.     emov (d1, v);
  600.       else
  601.     emov (d2, v);
  602.       break;
  603.     default:
  604.       emov (ezero, v);
  605.       break;
  606.     }
  607. PUT_REAL (v, value);
  608. }
  609.  
  610.  
  611. /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
  612.    implements REAL_VALUE_RNDZINT (x) (etrunci (x)).  */
  613.  
  614. REAL_VALUE_TYPE 
  615. etrunci (x)
  616.      REAL_VALUE_TYPE x;
  617. {
  618.   unsigned EMUSHORT f[NE], g[NE];
  619.   REAL_VALUE_TYPE r;
  620.   HOST_WIDE_INT l;
  621.  
  622.   GET_REAL (&x, g);
  623. #ifdef NANS
  624.   if (eisnan (g))
  625.     return (x);
  626. #endif
  627.   eifrac (g, &l, f);
  628.   ltoe (&l, g);
  629.   PUT_REAL (g, &r);
  630.   return (r);
  631. }
  632.  
  633.  
  634. /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
  635.    implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)).  */
  636.  
  637. REAL_VALUE_TYPE 
  638. etruncui (x)
  639.      REAL_VALUE_TYPE x;
  640. {
  641.   unsigned EMUSHORT f[NE], g[NE];
  642.   REAL_VALUE_TYPE r;
  643.   unsigned HOST_WIDE_INT l;
  644.  
  645.   GET_REAL (&x, g);
  646. #ifdef NANS
  647.   if (eisnan (g))
  648.     return (x);
  649. #endif
  650.   euifrac (g, &l, f);
  651.   ultoe (&l, g);
  652.   PUT_REAL (g, &r);
  653.   return (r);
  654. }
  655.  
  656.  
  657. /* This is the REAL_VALUE_ATOF function.  It converts a decimal string to
  658.    binary, rounding off as indicated by the machine_mode argument.  Then it
  659.    promotes the rounded value to REAL_VALUE_TYPE.  */
  660.  
  661. REAL_VALUE_TYPE 
  662. ereal_atof (s, t)
  663.      char *s;
  664.      enum machine_mode t;
  665. {
  666.   unsigned EMUSHORT tem[NE], e[NE];
  667.   REAL_VALUE_TYPE r;
  668.  
  669.   switch (t)
  670.     {
  671.     case HFmode:
  672.     case SFmode:
  673.       asctoe24 (s, tem);
  674.       e24toe (tem, e);
  675.       break;
  676.     case DFmode:
  677.       asctoe53 (s, tem);
  678.       e53toe (tem, e);
  679.       break;
  680.     case XFmode:
  681.       asctoe64 (s, tem);
  682.       e64toe (tem, e);
  683.       break;
  684.     case TFmode:
  685.       asctoe113 (s, tem);
  686.       e113toe (tem, e);
  687.       break;
  688.     default:
  689.       asctoe (s, e);
  690.     }
  691.   PUT_REAL (e, &r);
  692.   return (r);
  693. }
  694.  
  695.  
  696. /* Expansion of REAL_NEGATE.  */
  697.  
  698. REAL_VALUE_TYPE 
  699. ereal_negate (x)
  700.      REAL_VALUE_TYPE x;
  701. {
  702.   unsigned EMUSHORT e[NE];
  703.   REAL_VALUE_TYPE r;
  704.  
  705.   GET_REAL (&x, e);
  706.   eneg (e);
  707.   PUT_REAL (e, &r);
  708.   return (r);
  709. }
  710.  
  711.  
  712. /* Round real toward zero to HOST_WIDE_INT;
  713.    implements REAL_VALUE_FIX (x).  */
  714.  
  715. HOST_WIDE_INT
  716. efixi (x)
  717.      REAL_VALUE_TYPE x;
  718. {
  719.   unsigned EMUSHORT f[NE], g[NE];
  720.   HOST_WIDE_INT l;
  721.  
  722.   GET_REAL (&x, f);
  723. #ifdef NANS
  724.   if (eisnan (f))
  725.     {
  726.       warning ("conversion from NaN to int");
  727.       return (-1);
  728.     }
  729. #endif
  730.   eifrac (f, &l, g);
  731.   return l;
  732. }
  733.  
  734. /* Round real toward zero to unsigned HOST_WIDE_INT
  735.    implements  REAL_VALUE_UNSIGNED_FIX (x).
  736.    Negative input returns zero.  */
  737.  
  738. unsigned HOST_WIDE_INT
  739. efixui (x)
  740.      REAL_VALUE_TYPE x;
  741. {
  742.   unsigned EMUSHORT f[NE], g[NE];
  743.   unsigned HOST_WIDE_INT l;
  744.  
  745.   GET_REAL (&x, f);
  746. #ifdef NANS
  747.   if (eisnan (f))
  748.     {
  749.       warning ("conversion from NaN to unsigned int");
  750.       return (-1);
  751.     }
  752. #endif
  753.   euifrac (f, &l, g);
  754.   return l;
  755. }
  756.  
  757.  
  758. /* REAL_VALUE_FROM_INT macro.  */
  759.  
  760. void 
  761. ereal_from_int (d, i, j)
  762.      REAL_VALUE_TYPE *d;
  763.      HOST_WIDE_INT i, j;
  764. {
  765.   unsigned EMUSHORT df[NE], dg[NE];
  766.   HOST_WIDE_INT low, high;
  767.   int sign;
  768.  
  769.   sign = 0;
  770.   low = i;
  771.   if ((high = j) < 0)
  772.     {
  773.       sign = 1;
  774.       /* complement and add 1 */
  775.       high = ~high;
  776.       if (low)
  777.     low = -low;
  778.       else
  779.     high += 1;
  780.     }
  781.   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  782.   ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
  783.   emul (dg, df, dg);
  784.   ultoe ((unsigned HOST_WIDE_INT *) &low, df);
  785.   eadd (df, dg, dg);
  786.   if (sign)
  787.     eneg (dg);
  788.   PUT_REAL (dg, d);
  789. }
  790.  
  791.  
  792. /* REAL_VALUE_FROM_UNSIGNED_INT macro.   */
  793.  
  794. void 
  795. ereal_from_uint (d, i, j)
  796.      REAL_VALUE_TYPE *d;
  797.      unsigned HOST_WIDE_INT i, j;
  798. {
  799.   unsigned EMUSHORT df[NE], dg[NE];
  800.   unsigned HOST_WIDE_INT low, high;
  801.  
  802.   low = i;
  803.   high = j;
  804.   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  805.   ultoe (&high, dg);
  806.   emul (dg, df, dg);
  807.   ultoe (&low, df);
  808.   eadd (df, dg, dg);
  809.   PUT_REAL (dg, d);
  810. }
  811.  
  812.  
  813. /* REAL_VALUE_TO_INT macro.  */
  814.  
  815. void 
  816. ereal_to_int (low, high, rr)
  817.      HOST_WIDE_INT *low, *high;
  818.      REAL_VALUE_TYPE rr;
  819. {
  820.   unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
  821.   int s;
  822.  
  823.   GET_REAL (&rr, d);
  824. #ifdef NANS
  825.   if (eisnan (d))
  826.     {
  827.       warning ("conversion from NaN to int");
  828.       *low = -1;
  829.       *high = -1;
  830.       return;
  831.     }
  832. #endif
  833.   /* convert positive value */
  834.   s = 0;
  835.   if (eisneg (d))
  836.     {
  837.       eneg (d);
  838.       s = 1;
  839.     }
  840.   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  841.   ediv (df, d, dg);        /* dg = d / 2^32 is the high word */
  842.   euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
  843.   emul (df, dh, dg);        /* fractional part is the low word */
  844.   euifrac (dg, (unsigned HOST_WIDE_INT *)low, dh);
  845.   if (s)
  846.     {
  847.       /* complement and add 1 */
  848.       *high = ~(*high);
  849.       if (*low)
  850.     *low = -(*low);
  851.       else
  852.     *high += 1;
  853.     }
  854. }
  855.  
  856.  
  857. /* REAL_VALUE_LDEXP macro.  */
  858.  
  859. REAL_VALUE_TYPE
  860. ereal_ldexp (x, n)
  861.      REAL_VALUE_TYPE x;
  862.      int n;
  863. {
  864.   unsigned EMUSHORT e[NE], y[NE];
  865.   REAL_VALUE_TYPE r;
  866.  
  867.   GET_REAL (&x, e);
  868. #ifdef NANS
  869.   if (eisnan (e))
  870.     return (x);
  871. #endif
  872.   eldexp (e, n, y);
  873.   PUT_REAL (y, &r);
  874.   return (r);
  875. }
  876.  
  877. /* These routines are conditionally compiled because functions
  878.    of the same names may be defined in fold-const.c.  */
  879.  
  880. #ifdef REAL_ARITHMETIC
  881.  
  882. /* Check for infinity in a REAL_VALUE_TYPE. */
  883.  
  884. int
  885. target_isinf (x)
  886.      REAL_VALUE_TYPE x;
  887. {
  888.   unsigned EMUSHORT e[NE];
  889.  
  890. #ifdef INFINITY
  891.   GET_REAL (&x, e);
  892.   return (eisinf (e));
  893. #else
  894.   return 0;
  895. #endif
  896. }
  897.  
  898. /* Check whether a REAL_VALUE_TYPE item is a NaN. */
  899.  
  900. int
  901. target_isnan (x)
  902.      REAL_VALUE_TYPE x;
  903. {
  904.   unsigned EMUSHORT e[NE];
  905.  
  906. #ifdef NANS
  907.   GET_REAL (&x, e);
  908.   return (eisnan (e));
  909. #else
  910.   return (0);
  911. #endif
  912. }
  913.  
  914.  
  915. /* Check for a negative REAL_VALUE_TYPE number.
  916.    This just checks the sign bit, so that -0 counts as negative. */
  917.  
  918. int
  919. target_negative (x)
  920.      REAL_VALUE_TYPE x;
  921. {
  922.   return ereal_isneg (x);
  923. }
  924.  
  925. /* Expansion of REAL_VALUE_TRUNCATE.
  926.    The result is in floating point, rounded to nearest or even.  */
  927.  
  928. REAL_VALUE_TYPE
  929. real_value_truncate (mode, arg)
  930.      enum machine_mode mode;
  931.      REAL_VALUE_TYPE arg;
  932. {
  933.   unsigned EMUSHORT e[NE], t[NE];
  934.   REAL_VALUE_TYPE r;
  935.  
  936.   GET_REAL (&arg, e);
  937. #ifdef NANS
  938.   if (eisnan (e))
  939.     return (arg);
  940. #endif
  941.   eclear (t);
  942.   switch (mode)
  943.     {
  944.     case TFmode:
  945.       etoe113 (e, t);
  946.       e113toe (t, t);
  947.       break;
  948.  
  949.     case XFmode:
  950.       etoe64 (e, t);
  951.       e64toe (t, t);
  952.       break;
  953.  
  954.     case DFmode:
  955.       etoe53 (e, t);
  956.       e53toe (t, t);
  957.       break;
  958.  
  959.     case HFmode:
  960.     case SFmode:
  961.       etoe24 (e, t);
  962.       e24toe (t, t);
  963.       break;
  964.  
  965.     case SImode:
  966.       r = etrunci (arg);
  967.       return (r);
  968.  
  969.     /* If an unsupported type was requested, presume that
  970.        the machine files know something useful to do with
  971.        the unmodified value.  */
  972.  
  973.     default:
  974.       return (arg);
  975.     }
  976.   PUT_REAL (t, &r);
  977.   return (r);
  978. }
  979.  
  980. #endif /* REAL_ARITHMETIC defined */
  981.  
  982. /* Used for debugging--print the value of R in human-readable format
  983.    on stderr.  */
  984.  
  985. void
  986. debug_real (r)
  987.      REAL_VALUE_TYPE r;
  988. {
  989.   char dstr[30];
  990.  
  991.   REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
  992.   fprintf (stderr, "%s", dstr);
  993. }  
  994.  
  995.  
  996. /* The following routines convert REAL_VALUE_TYPE to the various floating
  997.    point formats that are meaningful to supported computers.
  998.  
  999.    The results are returned in 32-bit pieces, each piece stored in a `long'.  
  1000.    This is so they can be printed by statements like
  1001.  
  1002.       fprintf (file, "%lx, %lx", L[0],  L[1]);
  1003.  
  1004.    that will work on both narrow- and wide-word host computers.  */
  1005.  
  1006. /* Convert R to a 128-bit long double precision value.  The output array L
  1007.    contains four 32-bit pieces of the result, in the order they would appear
  1008.    in memory.  */
  1009.  
  1010. void 
  1011. etartdouble (r, l)
  1012.      REAL_VALUE_TYPE r;
  1013.      long l[];
  1014. {
  1015.   unsigned EMUSHORT e[NE];
  1016.  
  1017.   GET_REAL (&r, e);
  1018.   etoe113 (e, e);
  1019.   endian (e, l, TFmode);
  1020. }
  1021.  
  1022. /* Convert R to a double extended precision value.  The output array L
  1023.    contains three 32-bit pieces of the result, in the order they would
  1024.    appear in memory.  */
  1025.  
  1026. void 
  1027. etarldouble (r, l)
  1028.      REAL_VALUE_TYPE r;
  1029.      long l[];
  1030. {
  1031.   unsigned EMUSHORT e[NE];
  1032.  
  1033.   GET_REAL (&r, e);
  1034.   etoe64 (e, e);
  1035.   endian (e, l, XFmode);
  1036. }
  1037.  
  1038. /* Convert R to a double precision value.  The output array L contains two
  1039.    32-bit pieces of the result, in the order they would appear in memory.  */
  1040.  
  1041. void 
  1042. etardouble (r, l)
  1043.      REAL_VALUE_TYPE r;
  1044.      long l[];
  1045. {
  1046.   unsigned EMUSHORT e[NE];
  1047.  
  1048.   GET_REAL (&r, e);
  1049.   etoe53 (e, e);
  1050.   endian (e, l, DFmode);
  1051. }
  1052.  
  1053. /* Convert R to a single precision float value stored in the least-significant
  1054.    bits of a `long'.  */
  1055.  
  1056. long
  1057. etarsingle (r)
  1058.      REAL_VALUE_TYPE r;
  1059. {
  1060.   unsigned EMUSHORT e[NE];
  1061.   long l;
  1062.  
  1063.   GET_REAL (&r, e);
  1064.   etoe24 (e, e);
  1065.   endian (e, &l, SFmode);
  1066.   return ((long) l);
  1067. }
  1068.  
  1069. /* Convert X to a decimal ASCII string S for output to an assembly
  1070.    language file.  Note, there is no standard way to spell infinity or
  1071.    a NaN, so these values may require special treatment in the tm.h
  1072.    macros.  */
  1073.  
  1074. void
  1075. ereal_to_decimal (x, s)
  1076.      REAL_VALUE_TYPE x;
  1077.      char *s;
  1078. {
  1079.   unsigned EMUSHORT e[NE];
  1080.  
  1081.   GET_REAL (&x, e);
  1082.   etoasc (e, s, 20);
  1083. }
  1084.  
  1085. /* Compare X and Y.  Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
  1086.    or -2 if either is a NaN.   */
  1087.  
  1088. int
  1089. ereal_cmp (x, y)
  1090.      REAL_VALUE_TYPE x, y;
  1091. {
  1092.   unsigned EMUSHORT ex[NE], ey[NE];
  1093.  
  1094.   GET_REAL (&x, ex);
  1095.   GET_REAL (&y, ey);
  1096.   return (ecmp (ex, ey));
  1097. }
  1098.  
  1099. /*  Return 1 if the sign bit of X is set, else return 0.  */
  1100.  
  1101. int
  1102. ereal_isneg (x)
  1103.      REAL_VALUE_TYPE x;
  1104. {
  1105.   unsigned EMUSHORT ex[NE];
  1106.  
  1107.   GET_REAL (&x, ex);
  1108.   return (eisneg (ex));
  1109. }
  1110.  
  1111. /* End of REAL_ARITHMETIC interface */
  1112.  
  1113. /*
  1114.   Extended precision IEEE binary floating point arithmetic routines
  1115.  
  1116.   Numbers are stored in C language as arrays of 16-bit unsigned
  1117.   short integers.  The arguments of the routines are pointers to
  1118.   the arrays.
  1119.  
  1120.   External e type data structure, similar to Intel 8087 chip
  1121.   temporary real format but possibly with a larger significand:
  1122.  
  1123.     NE-1 significand words    (least significant word first,
  1124.                  most significant bit is normally set)
  1125.     exponent        (value = EXONE for 1.0,
  1126.                 top bit is the sign)
  1127.  
  1128.  
  1129.   Internal exploded e-type data structure of a number (a "word" is 16 bits):
  1130.  
  1131.   ei[0]    sign word    (0 for positive, 0xffff for negative)
  1132.   ei[1]    biased exponent    (value = EXONE for the number 1.0)
  1133.   ei[2]    high guard word    (always zero after normalization)
  1134.   ei[3]
  1135.   to ei[NI-2]    significand    (NI-4 significand words,
  1136.                   most significant word first,
  1137.                   most significant bit is set)
  1138.   ei[NI-1]    low guard word    (0x8000 bit is rounding place)
  1139.  
  1140.  
  1141.  
  1142.          Routines for external format e-type numbers
  1143.  
  1144.      asctoe (string, e)    ASCII string to extended double e type
  1145.      asctoe64 (string, &d)    ASCII string to long double
  1146.      asctoe53 (string, &d)    ASCII string to double
  1147.      asctoe24 (string, &f)    ASCII string to single
  1148.      asctoeg (string, e, prec) ASCII string to specified precision
  1149.      e24toe (&f, e)        IEEE single precision to e type
  1150.      e53toe (&d, e)        IEEE double precision to e type
  1151.      e64toe (&d, e)        IEEE long double precision to e type
  1152.      e113toe (&d, e)        128-bit long double precision to e type
  1153.      eabs (e)            absolute value
  1154.      eadd (a, b, c)        c = b + a
  1155.      eclear (e)        e = 0
  1156.      ecmp (a, b)        Returns 1 if a > b, 0 if a == b,
  1157.                  -1 if a < b, -2 if either a or b is a NaN.
  1158.      ediv (a, b, c)        c = b / a
  1159.      efloor (a, b)        truncate to integer, toward -infinity
  1160.      efrexp (a, exp, s)    extract exponent and significand
  1161.      eifrac (e, &l, frac)    e to HOST_WIDE_INT and e type fraction
  1162.      euifrac (e, &l, frac)   e to unsigned HOST_WIDE_INT and e type fraction
  1163.      einfin (e)        set e to infinity, leaving its sign alone
  1164.      eldexp (a, n, b)    multiply by 2**n
  1165.      emov (a, b)        b = a
  1166.      emul (a, b, c)        c = b * a
  1167.      eneg (e)            e = -e
  1168.      eround (a, b)        b = nearest integer value to a
  1169.      esub (a, b, c)        c = b - a
  1170.      e24toasc (&f, str, n)    single to ASCII string, n digits after decimal
  1171.      e53toasc (&d, str, n)    double to ASCII string, n digits after decimal
  1172.      e64toasc (&d, str, n)    80-bit long double to ASCII string
  1173.      e113toasc (&d, str, n)    128-bit long double to ASCII string
  1174.      etoasc (e, str, n)    e to ASCII string, n digits after decimal
  1175.      etoe24 (e, &f)        convert e type to IEEE single precision
  1176.      etoe53 (e, &d)        convert e type to IEEE double precision
  1177.      etoe64 (e, &d)        convert e type to IEEE long double precision
  1178.      ltoe (&l, e)        HOST_WIDE_INT to e type
  1179.      ultoe (&l, e)        unsigned HOST_WIDE_INT to e type
  1180.     eisneg (e)              1 if sign bit of e != 0, else 0
  1181.     eisinf (e)              1 if e has maximum exponent (non-IEEE)
  1182.                  or is infinite (IEEE)
  1183.         eisnan (e)              1 if e is a NaN
  1184.  
  1185.  
  1186.          Routines for internal format exploded e-type numbers
  1187.  
  1188.      eaddm (ai, bi)        add significands, bi = bi + ai
  1189.      ecleaz (ei)        ei = 0
  1190.      ecleazs (ei)        set ei = 0 but leave its sign alone
  1191.      ecmpm (ai, bi)        compare significands, return 1, 0, or -1
  1192.      edivm (ai, bi)        divide  significands, bi = bi / ai
  1193.      emdnorm (ai,l,s,exp)    normalize and round off
  1194.      emovi (a, ai)        convert external a to internal ai
  1195.      emovo (ai, a)        convert internal ai to external a
  1196.      emovz (ai, bi)        bi = ai, low guard word of bi = 0
  1197.      emulm (ai, bi)        multiply significands, bi = bi * ai
  1198.      enormlz (ei)        left-justify the significand
  1199.      eshdn1 (ai)        shift significand and guards down 1 bit
  1200.      eshdn8 (ai)        shift down 8 bits
  1201.      eshdn6 (ai)        shift down 16 bits
  1202.      eshift (ai, n)        shift ai n bits up (or down if n < 0)
  1203.      eshup1 (ai)        shift significand and guards up 1 bit
  1204.      eshup8 (ai)        shift up 8 bits
  1205.      eshup6 (ai)        shift up 16 bits
  1206.      esubm (ai, bi)        subtract significands, bi = bi - ai
  1207.         eiisinf (ai)            1 if infinite
  1208.         eiisnan (ai)            1 if a NaN
  1209.      eiisneg (ai)        1 if sign bit of ai != 0, else 0
  1210.         einan (ai)              set ai = NaN
  1211.         eiinfin (ai)            set ai = infinity
  1212.  
  1213.   The result is always normalized and rounded to NI-4 word precision
  1214.   after each arithmetic operation.
  1215.  
  1216.   Exception flags are NOT fully supported.
  1217.  
  1218.   Signaling NaN's are NOT supported; they are treated the same
  1219.   as quiet NaN's.
  1220.  
  1221.   Define INFINITY for support of infinity; otherwise a
  1222.   saturation arithmetic is implemented.
  1223.  
  1224.   Define NANS for support of Not-a-Number items; otherwise the
  1225.   arithmetic will never produce a NaN output, and might be confused
  1226.   by a NaN input.
  1227.   If NaN's are supported, the output of `ecmp (a,b)' is -2 if
  1228.   either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
  1229.   may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
  1230.   if in doubt.
  1231.  
  1232.   Denormals are always supported here where appropriate (e.g., not
  1233.   for conversion to DEC numbers).  */
  1234.  
  1235. /* Definitions for error codes that are passed to the common error handling
  1236.    routine mtherr.
  1237.  
  1238.    For Digital Equipment PDP-11 and VAX computers, certain
  1239.   IBM systems, and others that use numbers with a 56-bit
  1240.   significand, the symbol DEC should be defined.  In this
  1241.   mode, most floating point constants are given as arrays
  1242.   of octal integers to eliminate decimal to binary conversion
  1243.   errors that might be introduced by the compiler.
  1244.  
  1245.   For computers, such as IBM PC, that follow the IEEE
  1246.   Standard for Binary Floating Point Arithmetic (ANSI/IEEE
  1247.   Std 754-1985), the symbol IEEE should be defined.
  1248.   These numbers have 53-bit significands.  In this mode, constants
  1249.   are provided as arrays of hexadecimal 16 bit integers.
  1250.   The endian-ness of generated values is controlled by
  1251.   REAL_WORDS_BIG_ENDIAN.
  1252.  
  1253.   To accommodate other types of computer arithmetic, all
  1254.   constants are also provided in a normal decimal radix
  1255.   which one can hope are correctly converted to a suitable
  1256.   format by the available C language compiler.  To invoke
  1257.   this mode, the symbol UNK is defined.
  1258.  
  1259.   An important difference among these modes is a predefined
  1260.   set of machine arithmetic constants for each.  The numbers
  1261.   MACHEP (the machine roundoff error), MAXNUM (largest number
  1262.   represented), and several other parameters are preset by
  1263.   the configuration symbol.  Check the file const.c to
  1264.   ensure that these values are correct for your computer.
  1265.  
  1266.   For ANSI C compatibility, define ANSIC equal to 1.  Currently
  1267.   this affects only the atan2 function and others that use it. */
  1268.  
  1269. /* Constant definitions for math error conditions.  */
  1270.  
  1271. #define DOMAIN        1    /* argument domain error */
  1272. #define SING        2    /* argument singularity */
  1273. #define OVERFLOW    3    /* overflow range error */
  1274. #define UNDERFLOW    4    /* underflow range error */
  1275. #define TLOSS        5    /* total loss of precision */
  1276. #define PLOSS        6    /* partial loss of precision */
  1277. #define INVALID        7    /* NaN-producing operation */
  1278.  
  1279. /*  e type constants used by high precision check routines */
  1280.  
  1281. #if LONG_DOUBLE_TYPE_SIZE == 128
  1282. /* 0.0 */
  1283. unsigned EMUSHORT ezero[NE] =
  1284.  {0x0000, 0x0000, 0x0000, 0x0000,
  1285.   0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
  1286. extern unsigned EMUSHORT ezero[];
  1287.  
  1288. /* 5.0E-1 */
  1289. unsigned EMUSHORT ehalf[NE] =
  1290.  {0x0000, 0x0000, 0x0000, 0x0000,
  1291.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
  1292. extern unsigned EMUSHORT ehalf[];
  1293.  
  1294. /* 1.0E0 */
  1295. unsigned EMUSHORT eone[NE] =
  1296.  {0x0000, 0x0000, 0x0000, 0x0000,
  1297.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
  1298. extern unsigned EMUSHORT eone[];
  1299.  
  1300. /* 2.0E0 */
  1301. unsigned EMUSHORT etwo[NE] =
  1302.  {0x0000, 0x0000, 0x0000, 0x0000,
  1303.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
  1304. extern unsigned EMUSHORT etwo[];
  1305.  
  1306. /* 3.2E1 */
  1307. unsigned EMUSHORT e32[NE] =
  1308.  {0x0000, 0x0000, 0x0000, 0x0000,
  1309.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
  1310. extern unsigned EMUSHORT e32[];
  1311.  
  1312. /* 6.93147180559945309417232121458176568075500134360255E-1 */
  1313. unsigned EMUSHORT elog2[NE] =
  1314.  {0x40f3, 0xf6af, 0x03f2, 0xb398,
  1315.   0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
  1316. extern unsigned EMUSHORT elog2[];
  1317.  
  1318. /* 1.41421356237309504880168872420969807856967187537695E0 */
  1319. unsigned EMUSHORT esqrt2[NE] =
  1320.  {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
  1321.   0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
  1322. extern unsigned EMUSHORT esqrt2[];
  1323.  
  1324. /* 3.14159265358979323846264338327950288419716939937511E0 */
  1325. unsigned EMUSHORT epi[NE] =
  1326.  {0x2902, 0x1cd1, 0x80dc, 0x628b,
  1327.   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
  1328. extern unsigned EMUSHORT epi[];
  1329.  
  1330. #else
  1331. /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
  1332. unsigned EMUSHORT ezero[NE] =
  1333.  {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
  1334. unsigned EMUSHORT ehalf[NE] =
  1335.  {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
  1336. unsigned EMUSHORT eone[NE] =
  1337.  {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
  1338. unsigned EMUSHORT etwo[NE] =
  1339.  {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
  1340. unsigned EMUSHORT e32[NE] =
  1341.  {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
  1342. unsigned EMUSHORT elog2[NE] =
  1343.  {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
  1344. unsigned EMUSHORT esqrt2[NE] =
  1345.  {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
  1346. unsigned EMUSHORT epi[NE] =
  1347.  {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
  1348. #endif
  1349.  
  1350. /* Control register for rounding precision.
  1351.    This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits.  */
  1352.  
  1353. int rndprc = NBITS;
  1354. extern int rndprc;
  1355.  
  1356. /*  Clear out entire e-type number X.  */
  1357.  
  1358. static void 
  1359. eclear (x)
  1360.      register unsigned EMUSHORT *x;
  1361. {
  1362.   register int i;
  1363.  
  1364.   for (i = 0; i < NE; i++)
  1365.     *x++ = 0;
  1366. }
  1367.  
  1368. /* Move e-type number from A to B.  */
  1369.  
  1370. static void 
  1371. emov (a, b)
  1372.      register unsigned EMUSHORT *a, *b;
  1373. {
  1374.   register int i;
  1375.  
  1376.   for (i = 0; i < NE; i++)
  1377.     *b++ = *a++;
  1378. }
  1379.  
  1380.  
  1381. /* Absolute value of e-type X.  */
  1382.  
  1383. static void 
  1384. eabs (x)
  1385.      unsigned EMUSHORT x[];
  1386. {
  1387.   /* sign is top bit of last word of external format */
  1388.   x[NE - 1] &= 0x7fff;        
  1389. }
  1390.  
  1391. /* Negate the e-type number X.  */
  1392.  
  1393. static void 
  1394. eneg (x)
  1395.      unsigned EMUSHORT x[];
  1396. {
  1397.  
  1398.   x[NE - 1] ^= 0x8000;        /* Toggle the sign bit */
  1399. }
  1400.  
  1401. /* Return 1 if sign bit of e-type number X is nonzero, else zero.  */
  1402.  
  1403. static int 
  1404. eisneg (x)
  1405.      unsigned EMUSHORT x[];
  1406. {
  1407.  
  1408.   if (x[NE - 1] & 0x8000)
  1409.     return (1);
  1410.   else
  1411.     return (0);
  1412. }
  1413.  
  1414. /* Return 1 if e-type number X is infinity, else return zero.  */
  1415.  
  1416. static int 
  1417. eisinf (x)
  1418.      unsigned EMUSHORT x[];
  1419. {
  1420.  
  1421. #ifdef NANS
  1422.   if (eisnan (x))
  1423.     return (0);
  1424. #endif
  1425.   if ((x[NE - 1] & 0x7fff) == 0x7fff)
  1426.     return (1);
  1427.   else
  1428.     return (0);
  1429. }
  1430.  
  1431. /* Check if e-type number is not a number.  The bit pattern is one that we
  1432.    defined, so we know for sure how to detect it.  */
  1433.  
  1434. static int 
  1435. eisnan (x)
  1436.      unsigned EMUSHORT x[];
  1437. {
  1438. #ifdef NANS
  1439.   int i;
  1440.  
  1441.   /* NaN has maximum exponent */
  1442.   if ((x[NE - 1] & 0x7fff) != 0x7fff)
  1443.     return (0);
  1444.   /* ... and non-zero significand field. */
  1445.   for (i = 0; i < NE - 1; i++)
  1446.     {
  1447.       if (*x++ != 0)
  1448.         return (1);
  1449.     }
  1450. #endif
  1451.  
  1452.   return (0);
  1453. }
  1454.  
  1455. /*  Fill e-type number X with infinity pattern (IEEE)
  1456.     or largest possible number (non-IEEE). */
  1457.  
  1458. static void 
  1459. einfin (x)
  1460.      register unsigned EMUSHORT *x;
  1461. {
  1462.   register int i;
  1463.  
  1464. #ifdef INFINITY
  1465.   for (i = 0; i < NE - 1; i++)
  1466.     *x++ = 0;
  1467.   *x |= 32767;
  1468. #else
  1469.   for (i = 0; i < NE - 1; i++)
  1470.     *x++ = 0xffff;
  1471.   *x |= 32766;
  1472.   if (rndprc < NBITS)
  1473.     {
  1474.       if (rndprc == 113)
  1475.     {
  1476.       *(x - 9) = 0;
  1477.       *(x - 8) = 0;
  1478.     }
  1479.       if (rndprc == 64)
  1480.     {
  1481.       *(x - 5) = 0;
  1482.     }
  1483.       if (rndprc == 53)
  1484.     {
  1485.       *(x - 4) = 0xf800;
  1486.     }
  1487.       else
  1488.     {
  1489.       *(x - 4) = 0;
  1490.       *(x - 3) = 0;
  1491.       *(x - 2) = 0xff00;
  1492.     }
  1493.     }
  1494. #endif
  1495. }
  1496.  
  1497. /* Output an e-type NaN.
  1498.    This generates Intel's quiet NaN pattern for extended real.
  1499.    The exponent is 7fff, the leading mantissa word is c000.  */
  1500.  
  1501. static void 
  1502. enan (x, sign)
  1503.      register unsigned EMUSHORT *x;
  1504.      int sign;
  1505. {
  1506.   register int i;
  1507.  
  1508.   for (i = 0; i < NE - 2; i++)
  1509.     *x++ = 0;
  1510.   *x++ = 0xc000;
  1511.   *x = (sign << 15) | 0x7fff;
  1512. }
  1513.  
  1514. /* Move in an e-type number A, converting it to exploded e-type B.  */
  1515.  
  1516. static void 
  1517. emovi (a, b)
  1518.      unsigned EMUSHORT *a, *b;
  1519. {
  1520.   register unsigned EMUSHORT *p, *q;
  1521.   int i;
  1522.  
  1523.   q = b;
  1524.   p = a + (NE - 1);        /* point to last word of external number */
  1525.   /* get the sign bit */
  1526.   if (*p & 0x8000)
  1527.     *q++ = 0xffff;
  1528.   else
  1529.     *q++ = 0;
  1530.   /* get the exponent */
  1531.   *q = *p--;
  1532.   *q++ &= 0x7fff;        /* delete the sign bit */
  1533. #ifdef INFINITY
  1534.   if ((*(q - 1) & 0x7fff) == 0x7fff)
  1535.     {
  1536. #ifdef NANS
  1537.       if (eisnan (a))
  1538.     {
  1539.       *q++ = 0;
  1540.       for (i = 3; i < NI; i++)
  1541.         *q++ = *p--;
  1542.       return;
  1543.     }
  1544. #endif
  1545.  
  1546.       for (i = 2; i < NI; i++)
  1547.     *q++ = 0;
  1548.       return;
  1549.     }
  1550. #endif
  1551.  
  1552.   /* clear high guard word */
  1553.   *q++ = 0;
  1554.   /* move in the significand */
  1555.   for (i = 0; i < NE - 1; i++)
  1556.     *q++ = *p--;
  1557.   /* clear low guard word */
  1558.   *q = 0;
  1559. }
  1560.  
  1561. /* Move out exploded e-type number A, converting it to e type B.  */
  1562.  
  1563. static void 
  1564. emovo (a, b)
  1565.      unsigned EMUSHORT *a, *b;
  1566. {
  1567.   register unsigned EMUSHORT *p, *q;
  1568.   unsigned EMUSHORT i;
  1569.   int j;
  1570.  
  1571.   p = a;
  1572.   q = b + (NE - 1);        /* point to output exponent */
  1573.   /* combine sign and exponent */
  1574.   i = *p++;
  1575.   if (i)
  1576.     *q-- = *p++ | 0x8000;
  1577.   else
  1578.     *q-- = *p++;
  1579. #ifdef INFINITY
  1580.   if (*(p - 1) == 0x7fff)
  1581.     {
  1582. #ifdef NANS
  1583.       if (eiisnan (a))
  1584.     {
  1585.       enan (b, eiisneg (a));
  1586.       return;
  1587.     }
  1588. #endif
  1589.       einfin (b);
  1590.     return;
  1591.     }
  1592. #endif
  1593.   /* skip over guard word */
  1594.   ++p;
  1595.   /* move the significand */
  1596.   for (j = 0; j < NE - 1; j++)
  1597.     *q-- = *p++;
  1598. }
  1599.  
  1600. /* Clear out exploded e-type number XI.  */
  1601.  
  1602. static void 
  1603. ecleaz (xi)
  1604.      register unsigned EMUSHORT *xi;
  1605. {
  1606.   register int i;
  1607.  
  1608.   for (i = 0; i < NI; i++)
  1609.     *xi++ = 0;
  1610. }
  1611.  
  1612. /* Clear out exploded e-type XI, but don't touch the sign. */
  1613.  
  1614. static void 
  1615. ecleazs (xi)
  1616.      register unsigned EMUSHORT *xi;
  1617. {
  1618.   register int i;
  1619.  
  1620.   ++xi;
  1621.   for (i = 0; i < NI - 1; i++)
  1622.     *xi++ = 0;
  1623. }
  1624.  
  1625. /* Move exploded e-type number from A to B.  */
  1626.  
  1627. static void 
  1628. emovz (a, b)
  1629.      register unsigned EMUSHORT *a, *b;
  1630. {
  1631.   register int i;
  1632.  
  1633.   for (i = 0; i < NI - 1; i++)
  1634.     *b++ = *a++;
  1635.   /* clear low guard word */
  1636.   *b = 0;
  1637. }
  1638.  
  1639. /* Generate exploded e-type NaN.
  1640.    The explicit pattern for this is maximum exponent and
  1641.    top two significant bits set.  */
  1642.  
  1643. static void
  1644. einan (x)
  1645.      unsigned EMUSHORT x[];
  1646. {
  1647.  
  1648.   ecleaz (x);
  1649.   x[E] = 0x7fff;
  1650.   x[M + 1] = 0xc000;
  1651. }
  1652.  
  1653. /* Return nonzero if exploded e-type X is a NaN. */
  1654.  
  1655. static int 
  1656. eiisnan (x)
  1657.      unsigned EMUSHORT x[];
  1658. {
  1659.   int i;
  1660.  
  1661.   if ((x[E] & 0x7fff) == 0x7fff)
  1662.     {
  1663.       for (i = M + 1; i < NI; i++)
  1664.     {
  1665.       if (x[i] != 0)
  1666.         return (1);
  1667.     }
  1668.     }
  1669.   return (0);
  1670. }
  1671.  
  1672. /* Return nonzero if sign of exploded e-type X is nonzero.  */
  1673.  
  1674. static int 
  1675. eiisneg (x)
  1676.      unsigned EMUSHORT x[];
  1677. {
  1678.  
  1679.   return x[0] != 0;
  1680. }
  1681.  
  1682. /* Fill exploded e-type X with infinity pattern.
  1683.    This has maximum exponent and significand all zeros.  */
  1684.  
  1685. static void
  1686. eiinfin (x)
  1687.      unsigned EMUSHORT x[];
  1688. {
  1689.  
  1690.   ecleaz (x);
  1691.   x[E] = 0x7fff;
  1692. }
  1693.  
  1694. /* Return nonzero if exploded e-type X is infinite. */
  1695.  
  1696. static int 
  1697. eiisinf (x)
  1698.      unsigned EMUSHORT x[];
  1699. {
  1700.  
  1701. #ifdef NANS
  1702.   if (eiisnan (x))
  1703.     return (0);
  1704. #endif
  1705.   if ((x[E] & 0x7fff) == 0x7fff)
  1706.     return (1);
  1707.   return (0);
  1708. }
  1709.  
  1710.  
  1711. /* Compare significands of numbers in internal exploded e-type format.
  1712.    Guard words are included in the comparison.
  1713.  
  1714.    Returns    +1 if a > b
  1715.          0 if a == b
  1716.         -1 if a < b   */
  1717.  
  1718. static int
  1719. ecmpm (a, b)
  1720.      register unsigned EMUSHORT *a, *b;
  1721. {
  1722.   int i;
  1723.  
  1724.   a += M;            /* skip up to significand area */
  1725.   b += M;
  1726.   for (i = M; i < NI; i++)
  1727.     {
  1728.       if (*a++ != *b++)
  1729.     goto difrnt;
  1730.     }
  1731.   return (0);
  1732.  
  1733.  difrnt:
  1734.   if (*(--a) > *(--b))
  1735.     return (1);
  1736.   else
  1737.     return (-1);
  1738. }
  1739.  
  1740. /* Shift significand of exploded e-type X down by 1 bit.  */
  1741.  
  1742. static void 
  1743. eshdn1 (x)
  1744.      register unsigned EMUSHORT *x;
  1745. {
  1746.   register unsigned EMUSHORT bits;
  1747.   int i;
  1748.  
  1749.   x += M;            /* point to significand area */
  1750.  
  1751.   bits = 0;
  1752.   for (i = M; i < NI; i++)
  1753.     {
  1754.       if (*x & 1)
  1755.     bits |= 1;
  1756.       *x >>= 1;
  1757.       if (bits & 2)
  1758.     *x |= 0x8000;
  1759.       bits <<= 1;
  1760.       ++x;
  1761.     }
  1762. }
  1763.  
  1764. /* Shift significand of exploded e-type X up by 1 bit.  */
  1765.  
  1766. static void 
  1767. eshup1 (x)
  1768.      register unsigned EMUSHORT *x;
  1769. {
  1770.   register unsigned EMUSHORT bits;
  1771.   int i;
  1772.  
  1773.   x += NI - 1;
  1774.   bits = 0;
  1775.  
  1776.   for (i = M; i < NI; i++)
  1777.     {
  1778.       if (*x & 0x8000)
  1779.     bits |= 1;
  1780.       *x <<= 1;
  1781.       if (bits & 2)
  1782.     *x |= 1;
  1783.       bits <<= 1;
  1784.       --x;
  1785.     }
  1786. }
  1787.  
  1788.  
  1789. /* Shift significand of exploded e-type X down by 8 bits.  */
  1790.  
  1791. static void 
  1792. eshdn8 (x)
  1793.      register unsigned EMUSHORT *x;
  1794. {
  1795.   register unsigned EMUSHORT newbyt, oldbyt;
  1796.   int i;
  1797.  
  1798.   x += M;
  1799.   oldbyt = 0;
  1800.   for (i = M; i < NI; i++)
  1801.     {
  1802.       newbyt = *x << 8;
  1803.       *x >>= 8;
  1804.       *x |= oldbyt;
  1805.       oldbyt = newbyt;
  1806.       ++x;
  1807.     }
  1808. }
  1809.  
  1810. /* Shift significand of exploded e-type X up by 8 bits.  */
  1811.  
  1812. static void 
  1813. eshup8 (x)
  1814.      register unsigned EMUSHORT *x;
  1815. {
  1816.   int i;
  1817.   register unsigned EMUSHORT newbyt, oldbyt;
  1818.  
  1819.   x += NI - 1;
  1820.   oldbyt = 0;
  1821.  
  1822.   for (i = M; i < NI; i++)
  1823.     {
  1824.       newbyt = *x >> 8;
  1825.       *x <<= 8;
  1826.       *x |= oldbyt;
  1827.       oldbyt = newbyt;
  1828.       --x;
  1829.     }
  1830. }
  1831.  
  1832. /* Shift significand of exploded e-type X up by 16 bits.  */
  1833.  
  1834. static void 
  1835. eshup6 (x)
  1836.      register unsigned EMUSHORT *x;
  1837. {
  1838.   int i;
  1839.   register unsigned EMUSHORT *p;
  1840.  
  1841.   p = x + M;
  1842.   x += M + 1;
  1843.  
  1844.   for (i = M; i < NI - 1; i++)
  1845.     *p++ = *x++;
  1846.  
  1847.   *p = 0;
  1848. }
  1849.  
  1850. /* Shift significand of exploded e-type X down by 16 bits.  */
  1851.  
  1852. static void 
  1853. eshdn6 (x)
  1854.      register unsigned EMUSHORT *x;
  1855. {
  1856.   int i;
  1857.   register unsigned EMUSHORT *p;
  1858.  
  1859.   x += NI - 1;
  1860.   p = x + 1;
  1861.  
  1862.   for (i = M; i < NI - 1; i++)
  1863.     *(--p) = *(--x);
  1864.  
  1865.   *(--p) = 0;
  1866. }
  1867.  
  1868. /* Add significands of exploded e-type X and Y.  X + Y replaces Y.  */
  1869.  
  1870. static void 
  1871. eaddm (x, y)
  1872.      unsigned EMUSHORT *x, *y;
  1873. {
  1874.   register unsigned EMULONG a;
  1875.   int i;
  1876.   unsigned int carry;
  1877.  
  1878.   x += NI - 1;
  1879.   y += NI - 1;
  1880.   carry = 0;
  1881.   for (i = M; i < NI; i++)
  1882.     {
  1883.       a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
  1884.       if (a & 0x10000)
  1885.     carry = 1;
  1886.       else
  1887.     carry = 0;
  1888.       *y = (unsigned EMUSHORT) a;
  1889.       --x;
  1890.       --y;
  1891.     }
  1892. }
  1893.  
  1894. /* Subtract significands of exploded e-type X and Y.  Y - X replaces Y.  */
  1895.  
  1896. static void 
  1897. esubm (x, y)
  1898.      unsigned EMUSHORT *x, *y;
  1899. {
  1900.   unsigned EMULONG a;
  1901.   int i;
  1902.   unsigned int carry;
  1903.  
  1904.   x += NI - 1;
  1905.   y += NI - 1;
  1906.   carry = 0;
  1907.   for (i = M; i < NI; i++)
  1908.     {
  1909.       a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
  1910.       if (a & 0x10000)
  1911.     carry = 1;
  1912.       else
  1913.     carry = 0;
  1914.       *y = (unsigned EMUSHORT) a;
  1915.       --x;
  1916.       --y;
  1917.     }
  1918. }
  1919.  
  1920.  
  1921. static unsigned EMUSHORT equot[NI];
  1922.  
  1923.  
  1924. #if 0
  1925. /* Radix 2 shift-and-add versions of multiply and divide  */
  1926.  
  1927.  
  1928. /* Divide significands */
  1929.  
  1930. int 
  1931. edivm (den, num)
  1932.      unsigned EMUSHORT den[], num[];
  1933. {
  1934.   int i;
  1935.   register unsigned EMUSHORT *p, *q;
  1936.   unsigned EMUSHORT j;
  1937.  
  1938.   p = &equot[0];
  1939.   *p++ = num[0];
  1940.   *p++ = num[1];
  1941.  
  1942.   for (i = M; i < NI; i++)
  1943.     {
  1944.       *p++ = 0;
  1945.     }
  1946.  
  1947.   /* Use faster compare and subtraction if denominator has only 15 bits of
  1948.      significance.  */
  1949.  
  1950.   p = &den[M + 2];
  1951.   if (*p++ == 0)
  1952.     {
  1953.       for (i = M + 3; i < NI; i++)
  1954.     {
  1955.       if (*p++ != 0)
  1956.         goto fulldiv;
  1957.     }
  1958.       if ((den[M + 1] & 1) != 0)
  1959.     goto fulldiv;
  1960.       eshdn1 (num);
  1961.       eshdn1 (den);
  1962.  
  1963.       p = &den[M + 1];
  1964.       q = &num[M + 1];
  1965.  
  1966.       for (i = 0; i < NBITS + 2; i++)
  1967.     {
  1968.       if (*p <= *q)
  1969.         {
  1970.           *q -= *p;
  1971.           j = 1;
  1972.         }
  1973.       else
  1974.         {
  1975.           j = 0;
  1976.         }
  1977.       eshup1 (equot);
  1978.       equot[NI - 2] |= j;
  1979.       eshup1 (num);
  1980.     }
  1981.       goto divdon;
  1982.     }
  1983.  
  1984.   /* The number of quotient bits to calculate is NBITS + 1 scaling guard
  1985.      bit + 1 roundoff bit.  */
  1986.  
  1987.  fulldiv:
  1988.  
  1989.   p = &equot[NI - 2];
  1990.   for (i = 0; i < NBITS + 2; i++)
  1991.     {
  1992.       if (ecmpm (den, num) <= 0)
  1993.     {
  1994.       esubm (den, num);
  1995.       j = 1;        /* quotient bit = 1 */
  1996.     }
  1997.       else
  1998.     j = 0;
  1999.       eshup1 (equot);
  2000.       *p |= j;
  2001.       eshup1 (num);
  2002.     }
  2003.  
  2004.  divdon:
  2005.  
  2006.   eshdn1 (equot);
  2007.   eshdn1 (equot);
  2008.  
  2009.   /* test for nonzero remainder after roundoff bit */
  2010.   p = &num[M];
  2011.   j = 0;
  2012.   for (i = M; i < NI; i++)
  2013.     {
  2014.       j |= *p++;
  2015.     }
  2016.   if (j)
  2017.     j = 1;
  2018.  
  2019.  
  2020.   for (i = 0; i < NI; i++)
  2021.     num[i] = equot[i];
  2022.   return ((int) j);
  2023. }
  2024.  
  2025.  
  2026. /* Multiply significands */
  2027. int 
  2028. emulm (a, b)
  2029.      unsigned EMUSHORT a[], b[];
  2030. {
  2031.   unsigned EMUSHORT *p, *q;
  2032.   int i, j, k;
  2033.  
  2034.   equot[0] = b[0];
  2035.   equot[1] = b[1];
  2036.   for (i = M; i < NI; i++)
  2037.     equot[i] = 0;
  2038.  
  2039.   p = &a[NI - 2];
  2040.   k = NBITS;
  2041.   while (*p == 0)        /* significand is not supposed to be zero */
  2042.     {
  2043.       eshdn6 (a);
  2044.       k -= 16;
  2045.     }
  2046.   if ((*p & 0xff) == 0)
  2047.     {
  2048.       eshdn8 (a);
  2049.       k -= 8;
  2050.     }
  2051.  
  2052.   q = &equot[NI - 1];
  2053.   j = 0;
  2054.   for (i = 0; i < k; i++)
  2055.     {
  2056.       if (*p & 1)
  2057.     eaddm (b, equot);
  2058.       /* remember if there were any nonzero bits shifted out */
  2059.       if (*q & 1)
  2060.     j |= 1;
  2061.       eshdn1 (a);
  2062.       eshdn1 (equot);
  2063.     }
  2064.  
  2065.   for (i = 0; i < NI; i++)
  2066.     b[i] = equot[i];
  2067.  
  2068.   /* return flag for lost nonzero bits */
  2069.   return (j);
  2070. }
  2071.  
  2072. #else
  2073.  
  2074. /* Radix 65536 versions of multiply and divide.  */
  2075.  
  2076. /* Multiply significand of e-type number B
  2077.    by 16-bit quantity A, return e-type result to C. */
  2078.  
  2079. static void
  2080. m16m (a, b, c)
  2081.      unsigned int a;
  2082.      unsigned EMUSHORT b[], c[];
  2083. {
  2084.   register unsigned EMUSHORT *pp;
  2085.   register unsigned EMULONG carry;
  2086.   unsigned EMUSHORT *ps;
  2087.   unsigned EMUSHORT p[NI];
  2088.   unsigned EMULONG aa, m;
  2089.   int i;
  2090.  
  2091.   aa = a;
  2092.   pp = &p[NI-2];
  2093.   *pp++ = 0;
  2094.   *pp = 0;
  2095.   ps = &b[NI-1];
  2096.  
  2097.   for (i=M+1; i<NI; i++)
  2098.     {
  2099.       if (*ps == 0)
  2100.     {
  2101.       --ps;
  2102.       --pp;
  2103.       *(pp-1) = 0;
  2104.     }
  2105.       else
  2106.     {
  2107.       m = (unsigned EMULONG) aa * *ps--;
  2108.       carry = (m & 0xffff) + *pp;
  2109.       *pp-- = (unsigned EMUSHORT)carry;
  2110.       carry = (carry >> 16) + (m >> 16) + *pp;
  2111.       *pp = (unsigned EMUSHORT)carry;
  2112.       *(pp-1) = carry >> 16;
  2113.     }
  2114.     }
  2115.   for (i=M; i<NI; i++)
  2116.     c[i] = p[i];
  2117. }
  2118.  
  2119. /* Divide significands of exploded e-types NUM / DEN.  Neither the
  2120.    numerator NUM nor the denominator DEN is permitted to have its high guard
  2121.    word nonzero.  */
  2122.  
  2123. static int
  2124. edivm (den, num)
  2125.      unsigned EMUSHORT den[], num[];
  2126. {
  2127.   int i;
  2128.   register unsigned EMUSHORT *p;
  2129.   unsigned EMULONG tnum;
  2130.   unsigned EMUSHORT j, tdenm, tquot;
  2131.   unsigned EMUSHORT tprod[NI+1];
  2132.  
  2133.   p = &equot[0];
  2134.   *p++ = num[0];
  2135.   *p++ = num[1];
  2136.  
  2137.   for (i=M; i<NI; i++)
  2138.     {
  2139.       *p++ = 0;
  2140.     }
  2141.   eshdn1 (num);
  2142.   tdenm = den[M+1];
  2143.   for (i=M; i<NI; i++)
  2144.     {
  2145.       /* Find trial quotient digit (the radix is 65536). */
  2146.       tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
  2147.  
  2148.       /* Do not execute the divide instruction if it will overflow. */
  2149.       if ((tdenm * 0xffffL) < tnum)
  2150.     tquot = 0xffff;
  2151.       else
  2152.     tquot = tnum / tdenm;
  2153.       /* Multiply denominator by trial quotient digit. */
  2154.       m16m ((unsigned int)tquot, den, tprod);
  2155.       /* The quotient digit may have been overestimated. */
  2156.       if (ecmpm (tprod, num) > 0)
  2157.     {
  2158.       tquot -= 1;
  2159.       esubm (den, tprod);
  2160.       if (ecmpm (tprod, num) > 0)
  2161.         {
  2162.           tquot -= 1;
  2163.           esubm (den, tprod);
  2164.         }
  2165.     }
  2166.       esubm (tprod, num);
  2167.       equot[i] = tquot;
  2168.       eshup6(num);
  2169.     }
  2170.   /* test for nonzero remainder after roundoff bit */
  2171.   p = &num[M];
  2172.   j = 0;
  2173.   for (i=M; i<NI; i++)
  2174.     {
  2175.       j |= *p++;
  2176.     }
  2177.   if (j)
  2178.     j = 1;
  2179.  
  2180.   for (i=0; i<NI; i++)
  2181.     num[i] = equot[i];
  2182.  
  2183.   return ((int)j);
  2184. }
  2185.  
  2186. /* Multiply significands of exploded e-type A and B, result in B.  */
  2187.  
  2188. static int
  2189. emulm (a, b)
  2190.      unsigned EMUSHORT a[], b[];
  2191. {
  2192.   unsigned EMUSHORT *p, *q;
  2193.   unsigned EMUSHORT pprod[NI];
  2194.   unsigned EMUSHORT j;
  2195.   int i;
  2196.  
  2197.   equot[0] = b[0];
  2198.   equot[1] = b[1];
  2199.   for (i=M; i<NI; i++)
  2200.     equot[i] = 0;
  2201.  
  2202.   j = 0;
  2203.   p = &a[NI-1];
  2204.   q = &equot[NI-1];
  2205.   for (i=M+1; i<NI; i++)
  2206.     {
  2207.       if (*p == 0)
  2208.     {
  2209.       --p;
  2210.     }
  2211.       else
  2212.     {
  2213.       m16m ((unsigned int) *p--, b, pprod);
  2214.       eaddm(pprod, equot);
  2215.     }
  2216.       j |= *q;
  2217.       eshdn6(equot);
  2218.     }
  2219.  
  2220.   for (i=0; i<NI; i++)
  2221.     b[i] = equot[i];
  2222.  
  2223.   /* return flag for lost nonzero bits */
  2224.   return ((int)j);
  2225. }
  2226. #endif
  2227.  
  2228.  
  2229. /* Normalize and round off.
  2230.  
  2231.   The internal format number to be rounded is S.
  2232.   Input LOST is 0 if the value is exact.  This is the so-called sticky bit.
  2233.  
  2234.   Input SUBFLG indicates whether the number was obtained
  2235.   by a subtraction operation.  In that case if LOST is nonzero
  2236.   then the number is slightly smaller than indicated.
  2237.  
  2238.   Input EXP is the biased exponent, which may be negative.
  2239.   the exponent field of S is ignored but is replaced by
  2240.   EXP as adjusted by normalization and rounding.
  2241.  
  2242.   Input RCNTRL is the rounding control.  If it is nonzero, the
  2243.   returned value will be rounded to RNDPRC bits.
  2244.  
  2245.   For future reference:  In order for emdnorm to round off denormal
  2246.    significands at the right point, the input exponent must be
  2247.    adjusted to be the actual value it would have after conversion to
  2248.    the final floating point type.  This adjustment has been
  2249.    implemented for all type conversions (etoe53, etc.) and decimal
  2250.    conversions, but not for the arithmetic functions (eadd, etc.). 
  2251.    Data types having standard 15-bit exponents are not affected by
  2252.    this, but SFmode and DFmode are affected. For example, ediv with
  2253.    rndprc = 24 will not round correctly to 24-bit precision if the
  2254.    result is denormal.   */
  2255.  
  2256. static int rlast = -1;
  2257. static int rw = 0;
  2258. static unsigned EMUSHORT rmsk = 0;
  2259. static unsigned EMUSHORT rmbit = 0;
  2260. static unsigned EMUSHORT rebit = 0;
  2261. static int re = 0;
  2262. static unsigned EMUSHORT rbit[NI];
  2263.  
  2264. static void 
  2265. emdnorm (s, lost, subflg, exp, rcntrl)
  2266.      unsigned EMUSHORT s[];
  2267.      int lost;
  2268.      int subflg;
  2269.      EMULONG exp;
  2270.      int rcntrl;
  2271. {
  2272.   int i, j;
  2273.   unsigned EMUSHORT r;
  2274.  
  2275.   /* Normalize */
  2276.   j = enormlz (s);
  2277.  
  2278.   /* a blank significand could mean either zero or infinity. */
  2279. #ifndef INFINITY
  2280.   if (j > NBITS)
  2281.     {
  2282.       ecleazs (s);
  2283.       return;
  2284.     }
  2285. #endif
  2286.   exp -= j;
  2287. #ifndef INFINITY
  2288.   if (exp >= 32767L)
  2289.     goto overf;
  2290. #else
  2291.   if ((j > NBITS) && (exp < 32767))
  2292.     {
  2293.       ecleazs (s);
  2294.       return;
  2295.     }
  2296. #endif
  2297.   if (exp < 0L)
  2298.     {
  2299.       if (exp > (EMULONG) (-NBITS - 1))
  2300.     {
  2301.       j = (int) exp;
  2302.       i = eshift (s, j);
  2303.       if (i)
  2304.         lost = 1;
  2305.     }
  2306.       else
  2307.     {
  2308.       ecleazs (s);
  2309.       return;
  2310.     }
  2311.     }
  2312.   /* Round off, unless told not to by rcntrl. */
  2313.   if (rcntrl == 0)
  2314.     goto mdfin;
  2315.   /* Set up rounding parameters if the control register changed. */
  2316.   if (rndprc != rlast)
  2317.     {
  2318.       ecleaz (rbit);
  2319.       switch (rndprc)
  2320.     {
  2321.     default:
  2322.     case NBITS:
  2323.       rw = NI - 1;        /* low guard word */
  2324.       rmsk = 0xffff;
  2325.       rmbit = 0x8000;
  2326.       re = rw - 1;
  2327.       rebit = 1;
  2328.       break;
  2329.     case 113:
  2330.       rw = 10;
  2331.       rmsk = 0x7fff;
  2332.       rmbit = 0x4000;
  2333.       rebit = 0x8000;
  2334.       re = rw;
  2335.       break;
  2336.     case 64:
  2337.       rw = 7;
  2338.       rmsk = 0xffff;
  2339.       rmbit = 0x8000;
  2340.       re = rw - 1;
  2341.       rebit = 1;
  2342.       break;
  2343.       /* For DEC or IBM arithmetic */
  2344.     case 56:
  2345.       rw = 6;
  2346.       rmsk = 0xff;
  2347.       rmbit = 0x80;
  2348.       rebit = 0x100;
  2349.       re = rw;
  2350.       break;
  2351.     case 53:
  2352.       rw = 6;
  2353.       rmsk = 0x7ff;
  2354.       rmbit = 0x0400;
  2355.       rebit = 0x800;
  2356.       re = rw;
  2357.       break;
  2358.     case 24:
  2359.       rw = 4;
  2360.       rmsk = 0xff;
  2361.       rmbit = 0x80;
  2362.       rebit = 0x100;
  2363.       re = rw;
  2364.       break;
  2365.     }
  2366.       rbit[re] = rebit;
  2367.       rlast = rndprc;
  2368.     }
  2369.  
  2370.   /* Shift down 1 temporarily if the data structure has an implied
  2371.      most significant bit and the number is denormal.
  2372.      Intel long double denormals also lose one bit of precision.  */
  2373.   if ((exp <= 0) && (rndprc != NBITS)
  2374.       && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
  2375.     {
  2376.       lost |= s[NI - 1] & 1;
  2377.       eshdn1 (s);
  2378.     }
  2379.   /* Clear out all bits below the rounding bit,
  2380.      remembering in r if any were nonzero.  */
  2381.   r = s[rw] & rmsk;
  2382.   if (rndprc < NBITS)
  2383.     {
  2384.       i = rw + 1;
  2385.       while (i < NI)
  2386.     {
  2387.       if (s[i])
  2388.         r |= 1;
  2389.       s[i] = 0;
  2390.       ++i;
  2391.     }
  2392.     }
  2393.   s[rw] &= ~rmsk;
  2394.   if ((r & rmbit) != 0)
  2395.     {
  2396.       if (r == rmbit)
  2397.     {
  2398.       if (lost == 0)
  2399.         {            /* round to even */
  2400.           if ((s[re] & rebit) == 0)
  2401.         goto mddone;
  2402.         }
  2403.       else
  2404.         {
  2405.           if (subflg != 0)
  2406.         goto mddone;
  2407.         }
  2408.     }
  2409.       eaddm (rbit, s);
  2410.     }
  2411.  mddone:
  2412. /* Undo the temporary shift for denormal values. */
  2413.   if ((exp <= 0) && (rndprc != NBITS)
  2414.       && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
  2415.     {
  2416.       eshup1 (s);
  2417.     }
  2418.   if (s[2] != 0)
  2419.     {                /* overflow on roundoff */
  2420.       eshdn1 (s);
  2421.       exp += 1;
  2422.     }
  2423.  mdfin:
  2424.   s[NI - 1] = 0;
  2425.   if (exp >= 32767L)
  2426.     {
  2427. #ifndef INFINITY
  2428.     overf:
  2429. #endif
  2430. #ifdef INFINITY
  2431.       s[1] = 32767;
  2432.       for (i = 2; i < NI - 1; i++)
  2433.     s[i] = 0;
  2434.       if (extra_warnings)
  2435.     warning ("floating point overflow");
  2436. #else
  2437.       s[1] = 32766;
  2438.       s[2] = 0;
  2439.       for (i = M + 1; i < NI - 1; i++)
  2440.     s[i] = 0xffff;
  2441.       s[NI - 1] = 0;
  2442.       if ((rndprc < 64) || (rndprc == 113))
  2443.     {
  2444.       s[rw] &= ~rmsk;
  2445.       if (rndprc == 24)
  2446.         {
  2447.           s[5] = 0;
  2448.           s[6] = 0;
  2449.         }
  2450.     }
  2451. #endif
  2452.       return;
  2453.     }
  2454.   if (exp < 0)
  2455.     s[1] = 0;
  2456.   else
  2457.     s[1] = (unsigned EMUSHORT) exp;
  2458. }
  2459.  
  2460. /*  Subtract.  C = B - A, all e type numbers.  */
  2461.  
  2462. static int subflg = 0;
  2463.  
  2464. static void 
  2465. esub (a, b, c)
  2466.      unsigned EMUSHORT *a, *b, *c;
  2467. {
  2468.  
  2469. #ifdef NANS
  2470.   if (eisnan (a))
  2471.     {
  2472.       emov (a, c);
  2473.       return;
  2474.     }
  2475.   if (eisnan (b))
  2476.     {
  2477.       emov (b, c);
  2478.       return;
  2479.     }
  2480. /* Infinity minus infinity is a NaN.
  2481.    Test for subtracting infinities of the same sign. */
  2482.   if (eisinf (a) && eisinf (b)
  2483.       && ((eisneg (a) ^ eisneg (b)) == 0))
  2484.     {
  2485.       mtherr ("esub", INVALID);
  2486.       enan (c, 0);
  2487.       return;
  2488.     }
  2489. #endif
  2490.   subflg = 1;
  2491.   eadd1 (a, b, c);
  2492. }
  2493.  
  2494. /* Add.  C = A + B, all e type. */
  2495.  
  2496. static void 
  2497. eadd (a, b, c)
  2498.      unsigned EMUSHORT *a, *b, *c;
  2499. {
  2500.  
  2501. #ifdef NANS
  2502. /* NaN plus anything is a NaN. */
  2503.   if (eisnan (a))
  2504.     {
  2505.       emov (a, c);
  2506.       return;
  2507.     }
  2508.   if (eisnan (b))
  2509.     {
  2510.       emov (b, c);
  2511.       return;
  2512.     }
  2513. /* Infinity minus infinity is a NaN.
  2514.    Test for adding infinities of opposite signs. */
  2515.   if (eisinf (a) && eisinf (b)
  2516.       && ((eisneg (a) ^ eisneg (b)) != 0))
  2517.     {
  2518.       mtherr ("esub", INVALID);
  2519.       enan (c, 0);
  2520.       return;
  2521.     }
  2522. #endif
  2523.   subflg = 0;
  2524.   eadd1 (a, b, c);
  2525. }
  2526.  
  2527. /* Arithmetic common to both addition and subtraction.  */
  2528.  
  2529. static void 
  2530. eadd1 (a, b, c)
  2531.      unsigned EMUSHORT *a, *b, *c;
  2532. {
  2533.   unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
  2534.   int i, lost, j, k;
  2535.   EMULONG lt, lta, ltb;
  2536.  
  2537. #ifdef INFINITY
  2538.   if (eisinf (a))
  2539.     {
  2540.       emov (a, c);
  2541.       if (subflg)
  2542.     eneg (c);
  2543.       return;
  2544.     }
  2545.   if (eisinf (b))
  2546.     {
  2547.       emov (b, c);
  2548.       return;
  2549.     }
  2550. #endif
  2551.   emovi (a, ai);
  2552.   emovi (b, bi);
  2553.   if (subflg)
  2554.     ai[0] = ~ai[0];
  2555.  
  2556.   /* compare exponents */
  2557.   lta = ai[E];
  2558.   ltb = bi[E];
  2559.   lt = lta - ltb;
  2560.   if (lt > 0L)
  2561.     {                /* put the larger number in bi */
  2562.       emovz (bi, ci);
  2563.       emovz (ai, bi);
  2564.       emovz (ci, ai);
  2565.       ltb = bi[E];
  2566.       lt = -lt;
  2567.     }
  2568.   lost = 0;
  2569.   if (lt != 0L)
  2570.     {
  2571.       if (lt < (EMULONG) (-NBITS - 1))
  2572.     goto done;        /* answer same as larger addend */
  2573.       k = (int) lt;
  2574.       lost = eshift (ai, k);    /* shift the smaller number down */
  2575.     }
  2576.   else
  2577.     {
  2578.       /* exponents were the same, so must compare significands */
  2579.       i = ecmpm (ai, bi);
  2580.       if (i == 0)
  2581.     {            /* the numbers are identical in magnitude */
  2582.       /* if different signs, result is zero */
  2583.       if (ai[0] != bi[0])
  2584.         {
  2585.           eclear (c);
  2586.           return;
  2587.         }
  2588.       /* if same sign, result is double */
  2589.       /* double denormalized tiny number */
  2590.       if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
  2591.         {
  2592.           eshup1 (bi);
  2593.           goto done;
  2594.         }
  2595.       /* add 1 to exponent unless both are zero! */
  2596.       for (j = 1; j < NI - 1; j++)
  2597.         {
  2598.           if (bi[j] != 0)
  2599.         {
  2600.           /* This could overflow, but let emovo take care of that. */
  2601.           ltb += 1;
  2602.           break;
  2603.         }
  2604.         }
  2605.       bi[E] = (unsigned EMUSHORT) ltb;
  2606.       goto done;
  2607.     }
  2608.       if (i > 0)
  2609.     {            /* put the larger number in bi */
  2610.       emovz (bi, ci);
  2611.       emovz (ai, bi);
  2612.       emovz (ci, ai);
  2613.     }
  2614.     }
  2615.   if (ai[0] == bi[0])
  2616.     {
  2617.       eaddm (ai, bi);
  2618.       subflg = 0;
  2619.     }
  2620.   else
  2621.     {
  2622.       esubm (ai, bi);
  2623.       subflg = 1;
  2624.     }
  2625.   emdnorm (bi, lost, subflg, ltb, 64);
  2626.  
  2627.  done:
  2628.   emovo (bi, c);
  2629. }
  2630.  
  2631. /* Divide: C = B/A, all e type.  */
  2632.  
  2633. static void 
  2634. ediv (a, b, c)
  2635.      unsigned EMUSHORT *a, *b, *c;
  2636. {
  2637.   unsigned EMUSHORT ai[NI], bi[NI];
  2638.   int i;
  2639.   EMULONG lt, lta, ltb;
  2640.  
  2641. #ifdef NANS
  2642. /* Return any NaN input. */
  2643.   if (eisnan (a))
  2644.     {
  2645.     emov (a, c);
  2646.     return;
  2647.     }
  2648.   if (eisnan (b))
  2649.     {
  2650.     emov (b, c);
  2651.     return;
  2652.     }
  2653. /* Zero over zero, or infinity over infinity, is a NaN. */
  2654.   if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
  2655.       || (eisinf (a) && eisinf (b)))
  2656.     {
  2657.     mtherr ("ediv", INVALID);
  2658.     enan (c, eisneg (a) ^ eisneg (b));
  2659.     return;
  2660.     }
  2661. #endif
  2662. /* Infinity over anything else is infinity. */
  2663. #ifdef INFINITY
  2664.   if (eisinf (b))
  2665.     {
  2666.       if (eisneg (a) ^ eisneg (b))
  2667.     *(c + (NE - 1)) = 0x8000;
  2668.       else
  2669.     *(c + (NE - 1)) = 0;
  2670.       einfin (c);
  2671.       return;
  2672.     }
  2673. /* Anything else over infinity is zero. */
  2674.   if (eisinf (a))
  2675.     {
  2676.       eclear (c);
  2677.       return;
  2678.     }
  2679. #endif
  2680.   emovi (a, ai);
  2681.   emovi (b, bi);
  2682.   lta = ai[E];
  2683.   ltb = bi[E];
  2684.   if (bi[E] == 0)
  2685.     {                /* See if numerator is zero. */
  2686.       for (i = 1; i < NI - 1; i++)
  2687.     {
  2688.       if (bi[i] != 0)
  2689.         {
  2690.           ltb -= enormlz (bi);
  2691.           goto dnzro1;
  2692.         }
  2693.     }
  2694.       eclear (c);
  2695.       return;
  2696.     }
  2697.  dnzro1:
  2698.  
  2699.   if (ai[E] == 0)
  2700.     {                /* possible divide by zero */
  2701.       for (i = 1; i < NI - 1; i++)
  2702.     {
  2703.       if (ai[i] != 0)
  2704.         {
  2705.           lta -= enormlz (ai);
  2706.           goto dnzro2;
  2707.         }
  2708.     }
  2709.       if (ai[0] == bi[0])
  2710.     *(c + (NE - 1)) = 0;
  2711.       else
  2712.     *(c + (NE - 1)) = 0x8000;
  2713. /* Divide by zero is not an invalid operation.
  2714.    It is a divide-by-zero operation!   */
  2715.       einfin (c);
  2716.       mtherr ("ediv", SING);
  2717.       return;
  2718.     }
  2719.  dnzro2:
  2720.  
  2721.   i = edivm (ai, bi);
  2722.   /* calculate exponent */
  2723.   lt = ltb - lta + EXONE;
  2724.   emdnorm (bi, i, 0, lt, 64);
  2725.   /* set the sign */
  2726.   if (ai[0] == bi[0])
  2727.     bi[0] = 0;
  2728.   else
  2729.     bi[0] = 0Xffff;
  2730.   emovo (bi, c);
  2731. }
  2732.  
  2733. /* Multiply e-types A and B, return e-type product C.   */
  2734.  
  2735. static void 
  2736. emul (a, b, c)
  2737.      unsigned EMUSHORT *a, *b, *c;
  2738. {
  2739.   unsigned EMUSHORT ai[NI], bi[NI];
  2740.   int i, j;
  2741.   EMULONG lt, lta, ltb;
  2742.  
  2743. #ifdef NANS
  2744. /* NaN times anything is the same NaN. */
  2745.   if (eisnan (a))
  2746.     {
  2747.     emov (a, c);
  2748.     return;
  2749.     }
  2750.   if (eisnan (b))
  2751.     {
  2752.     emov (b, c);
  2753.     return;
  2754.     }
  2755. /* Zero times infinity is a NaN. */
  2756.   if ((eisinf (a) && (ecmp (b, ezero) == 0))
  2757.       || (eisinf (b) && (ecmp (a, ezero) == 0)))
  2758.     {
  2759.     mtherr ("emul", INVALID);
  2760.     enan (c, eisneg (a) ^ eisneg (b));
  2761.     return;
  2762.     }
  2763. #endif
  2764. /* Infinity times anything else is infinity. */
  2765. #ifdef INFINITY
  2766.   if (eisinf (a) || eisinf (b))
  2767.     {
  2768.       if (eisneg (a) ^ eisneg (b))
  2769.     *(c + (NE - 1)) = 0x8000;
  2770.       else
  2771.     *(c + (NE - 1)) = 0;
  2772.       einfin (c);
  2773.       return;
  2774.     }
  2775. #endif
  2776.   emovi (a, ai);
  2777.   emovi (b, bi);
  2778.   lta = ai[E];
  2779.   ltb = bi[E];
  2780.   if (ai[E] == 0)
  2781.     {
  2782.       for (i = 1; i < NI - 1; i++)
  2783.     {
  2784.       if (ai[i] != 0)
  2785.         {
  2786.           lta -= enormlz (ai);
  2787.           goto mnzer1;
  2788.         }
  2789.     }
  2790.       eclear (c);
  2791.       return;
  2792.     }
  2793.  mnzer1:
  2794.  
  2795.   if (bi[E] == 0)
  2796.     {
  2797.       for (i = 1; i < NI - 1; i++)
  2798.     {
  2799.       if (bi[i] != 0)
  2800.         {
  2801.           ltb -= enormlz (bi);
  2802.           goto mnzer2;
  2803.         }
  2804.     }
  2805.       eclear (c);
  2806.       return;
  2807.     }
  2808.  mnzer2:
  2809.  
  2810.   /* Multiply significands */
  2811.   j = emulm (ai, bi);
  2812.   /* calculate exponent */
  2813.   lt = lta + ltb - (EXONE - 1);
  2814.   emdnorm (bi, j, 0, lt, 64);
  2815.   /* calculate sign of product */
  2816.   if (ai[0] == bi[0])
  2817.     bi[0] = 0;
  2818.   else
  2819.     bi[0] = 0xffff;
  2820.   emovo (bi, c);
  2821. }
  2822.  
  2823. /* Convert double precision PE to e-type Y.  */
  2824.  
  2825. static void
  2826. e53toe (pe, y)
  2827.      unsigned EMUSHORT *pe, *y;
  2828. {
  2829. #ifdef DEC
  2830.  
  2831.   dectoe (pe, y);
  2832.  
  2833. #else
  2834. #ifdef IBM
  2835.  
  2836.   ibmtoe (pe, y, DFmode);
  2837.  
  2838. #else
  2839.   register unsigned EMUSHORT r;
  2840.   register unsigned EMUSHORT *e, *p;
  2841.   unsigned EMUSHORT yy[NI];
  2842.   int denorm, k;
  2843.  
  2844.   e = pe;
  2845.   denorm = 0;            /* flag if denormalized number */
  2846.   ecleaz (yy);
  2847.   if (! REAL_WORDS_BIG_ENDIAN)
  2848.     e += 3;
  2849.   r = *e;
  2850.   yy[0] = 0;
  2851.   if (r & 0x8000)
  2852.     yy[0] = 0xffff;
  2853.   yy[M] = (r & 0x0f) | 0x10;
  2854.   r &= ~0x800f;            /* strip sign and 4 significand bits */
  2855. #ifdef INFINITY
  2856.   if (r == 0x7ff0)
  2857.     {
  2858. #ifdef NANS
  2859.       if (! REAL_WORDS_BIG_ENDIAN)
  2860.     {
  2861.       if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
  2862.           || (pe[1] != 0) || (pe[0] != 0))
  2863.         {
  2864.           enan (y, yy[0] != 0);
  2865.           return;
  2866.         }
  2867.     }
  2868.       else
  2869.     {
  2870.       if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
  2871.           || (pe[2] != 0) || (pe[3] != 0))
  2872.         {
  2873.           enan (y, yy[0] != 0);
  2874.           return;
  2875.         }
  2876.     }
  2877. #endif  /* NANS */
  2878.       eclear (y);
  2879.       einfin (y);
  2880.       if (yy[0])
  2881.     eneg (y);
  2882.       return;
  2883.     }
  2884. #endif  /* INFINITY */
  2885.   r >>= 4;
  2886.   /* If zero exponent, then the significand is denormalized.
  2887.      So take back the understood high significand bit. */
  2888.  
  2889.   if (r == 0)
  2890.     {
  2891.       denorm = 1;
  2892.       yy[M] &= ~0x10;
  2893.     }
  2894.   r += EXONE - 01777;
  2895.   yy[E] = r;
  2896.   p = &yy[M + 1];
  2897. #ifdef IEEE
  2898.   if (! REAL_WORDS_BIG_ENDIAN)
  2899.     {
  2900.       *p++ = *(--e);
  2901.       *p++ = *(--e);
  2902.       *p++ = *(--e);
  2903.     }
  2904.   else
  2905.     {
  2906.       ++e;
  2907.       *p++ = *e++;
  2908.       *p++ = *e++;
  2909.       *p++ = *e++;
  2910.     }
  2911. #endif
  2912.   eshift (yy, -5);
  2913.   if (denorm)
  2914.     {                /* if zero exponent, then normalize the significand */
  2915.       if ((k = enormlz (yy)) > NBITS)
  2916.     ecleazs (yy);
  2917.       else
  2918.     yy[E] -= (unsigned EMUSHORT) (k - 1);
  2919.     }
  2920.   emovo (yy, y);
  2921. #endif /* not IBM */
  2922. #endif /* not DEC */
  2923. }
  2924.  
  2925. /* Convert double extended precision float PE to e type Y.  */
  2926.  
  2927. static void 
  2928. e64toe (pe, y)
  2929.      unsigned EMUSHORT *pe, *y;
  2930. {
  2931.   unsigned EMUSHORT yy[NI];
  2932.   unsigned EMUSHORT *e, *p, *q;
  2933.   int i;
  2934.  
  2935.   e = pe;
  2936.   p = yy;
  2937.   for (i = 0; i < NE - 5; i++)
  2938.     *p++ = 0;
  2939. /* This precision is not ordinarily supported on DEC or IBM. */
  2940. #ifdef DEC
  2941.   for (i = 0; i < 5; i++)
  2942.     *p++ = *e++;
  2943. #endif
  2944. #ifdef IBM
  2945.   p = &yy[0] + (NE - 1);
  2946.   *p-- = *e++;
  2947.   ++e;
  2948.   for (i = 0; i < 5; i++)
  2949.     *p-- = *e++;
  2950. #endif
  2951. #ifdef IEEE
  2952.   if (! REAL_WORDS_BIG_ENDIAN)
  2953.     {
  2954.       for (i = 0; i < 5; i++)
  2955.     *p++ = *e++;
  2956.  
  2957.       /* For denormal long double Intel format, shift significand up one
  2958.      -- but only if the top significand bit is zero.  A top bit of 1
  2959.      is "pseudodenormal" when the exponent is zero.  */
  2960.       if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
  2961.     {
  2962.       unsigned EMUSHORT temp[NI];
  2963.  
  2964.       emovi(yy, temp);
  2965.       eshup1(temp);
  2966.       emovo(temp,y);
  2967.       return;
  2968.     }
  2969.     }
  2970.   else
  2971.     {
  2972.       p = &yy[0] + (NE - 1);
  2973.       *p-- = *e++;
  2974.       ++e;
  2975.       for (i = 0; i < 4; i++)
  2976.     *p-- = *e++;
  2977.     }
  2978. #endif
  2979. #ifdef INFINITY
  2980.   /* Point to the exponent field and check max exponent cases.  */
  2981.   p = &yy[NE - 1];
  2982.   if (*p == 0x7fff)
  2983.     {
  2984. #ifdef NANS
  2985.       if (! REAL_WORDS_BIG_ENDIAN)
  2986.     {
  2987.       for (i = 0; i < 4; i++)
  2988.         {
  2989.           if ((i != 3 && pe[i] != 0)
  2990.           /* Anything but 0x8000 here, including 0, is a NaN.  */
  2991.           || (i == 3 && pe[i] != 0x8000))
  2992.         {
  2993.           enan (y, (*p & 0x8000) != 0);
  2994.           return;
  2995.         }
  2996.         }
  2997.     }
  2998.       else
  2999.     {
  3000.       for (i = 1; i <= 4; i++)
  3001.         {
  3002.           if (pe[i] != 0)
  3003.         {
  3004.           enan (y, (*p & 0x8000) != 0);
  3005.           return;
  3006.         }
  3007.         }
  3008.     }
  3009. #endif /* NANS */
  3010.       eclear (y);
  3011.       einfin (y);
  3012.       if (*p & 0x8000)
  3013.     eneg (y);
  3014.       return;
  3015.     }
  3016. #endif  /* INFINITY */
  3017.   p = yy;
  3018.   q = y;
  3019.   for (i = 0; i < NE; i++)
  3020.     *q++ = *p++;
  3021. }
  3022.  
  3023. /* Convert 128-bit long double precision float PE to e type Y.  */
  3024.  
  3025. static void 
  3026. e113toe (pe, y)
  3027.      unsigned EMUSHORT *pe, *y;
  3028. {
  3029.   register unsigned EMUSHORT r;
  3030.   unsigned EMUSHORT *e, *p;
  3031.   unsigned EMUSHORT yy[NI];
  3032.   int denorm, i;
  3033.  
  3034.   e = pe;
  3035.   denorm = 0;
  3036.   ecleaz (yy);
  3037. #ifdef IEEE
  3038.   if (! REAL_WORDS_BIG_ENDIAN)
  3039.     e += 7;
  3040. #endif
  3041.   r = *e;
  3042.   yy[0] = 0;
  3043.   if (r & 0x8000)
  3044.     yy[0] = 0xffff;
  3045.   r &= 0x7fff;
  3046. #ifdef INFINITY
  3047.   if (r == 0x7fff)
  3048.     {
  3049. #ifdef NANS
  3050.       if (! REAL_WORDS_BIG_ENDIAN)
  3051.     {
  3052.       for (i = 0; i < 7; i++)
  3053.         {
  3054.           if (pe[i] != 0)
  3055.         {
  3056.           enan (y, yy[0] != 0);
  3057.           return;
  3058.         }
  3059.         }
  3060.     }
  3061.       else
  3062.     {
  3063.       for (i = 1; i < 8; i++)
  3064.         {
  3065.           if (pe[i] != 0)
  3066.         {
  3067.           enan (y, yy[0] != 0);
  3068.           return;
  3069.         }
  3070.         }
  3071.     }
  3072. #endif /* NANS */
  3073.       eclear (y);
  3074.       einfin (y);
  3075.       if (yy[0])
  3076.     eneg (y);
  3077.       return;
  3078.     }
  3079. #endif  /* INFINITY */
  3080.   yy[E] = r;
  3081.   p = &yy[M + 1];
  3082. #ifdef IEEE
  3083.   if (! REAL_WORDS_BIG_ENDIAN)
  3084.     {
  3085.       for (i = 0; i < 7; i++)
  3086.     *p++ = *(--e);
  3087.     }
  3088.   else
  3089.     {
  3090.       ++e;
  3091.       for (i = 0; i < 7; i++)
  3092.     *p++ = *e++;
  3093.     }
  3094. #endif
  3095. /* If denormal, remove the implied bit; else shift down 1. */
  3096.   if (r == 0)
  3097.     {
  3098.       yy[M] = 0;
  3099.     }
  3100.   else
  3101.     {
  3102.       yy[M] = 1;
  3103.       eshift (yy, -1);
  3104.     }
  3105.   emovo (yy, y);
  3106. }
  3107.  
  3108. /* Convert single precision float PE to e type Y.  */
  3109.  
  3110. static void 
  3111. e24toe (pe, y)
  3112.      unsigned EMUSHORT *pe, *y;
  3113. {
  3114. #ifdef IBM
  3115.  
  3116.   ibmtoe (pe, y, SFmode);
  3117.  
  3118. #else
  3119.   register unsigned EMUSHORT r;
  3120.   register unsigned EMUSHORT *e, *p;
  3121.   unsigned EMUSHORT yy[NI];
  3122.   int denorm, k;
  3123.  
  3124.   e = pe;
  3125.   denorm = 0;            /* flag if denormalized number */
  3126.   ecleaz (yy);
  3127. #ifdef IEEE
  3128.   if (! REAL_WORDS_BIG_ENDIAN)
  3129.     e += 1;
  3130. #endif
  3131. #ifdef DEC
  3132.   e += 1;
  3133. #endif
  3134.   r = *e;
  3135.   yy[0] = 0;
  3136.   if (r & 0x8000)
  3137.     yy[0] = 0xffff;
  3138.   yy[M] = (r & 0x7f) | 0200;
  3139.   r &= ~0x807f;            /* strip sign and 7 significand bits */
  3140. #ifdef INFINITY
  3141.   if (r == 0x7f80)
  3142.     {
  3143. #ifdef NANS
  3144.       if (REAL_WORDS_BIG_ENDIAN)
  3145.     {
  3146.       if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
  3147.         {
  3148.           enan (y, yy[0] != 0);
  3149.           return;
  3150.         }
  3151.     }
  3152.       else
  3153.     {
  3154.       if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
  3155.         {
  3156.           enan (y, yy[0] != 0);
  3157.           return;
  3158.         }
  3159.     }
  3160. #endif  /* NANS */
  3161.       eclear (y);
  3162.       einfin (y);
  3163.       if (yy[0])
  3164.     eneg (y);
  3165.       return;
  3166.     }
  3167. #endif  /* INFINITY */
  3168.   r >>= 7;
  3169.   /* If zero exponent, then the significand is denormalized.
  3170.      So take back the understood high significand bit. */
  3171.   if (r == 0)
  3172.     {
  3173.       denorm = 1;
  3174.       yy[M] &= ~0200;
  3175.     }
  3176.   r += EXONE - 0177;
  3177.   yy[E] = r;
  3178.   p = &yy[M + 1];
  3179. #ifdef DEC
  3180.   *p++ = *(--e);
  3181. #endif
  3182. #ifdef IEEE
  3183.   if (! REAL_WORDS_BIG_ENDIAN)
  3184.     *p++ = *(--e);
  3185.   else
  3186.     {
  3187.       ++e;
  3188.       *p++ = *e++;
  3189.     }
  3190. #endif
  3191.   eshift (yy, -8);
  3192.   if (denorm)
  3193.     {                /* if zero exponent, then normalize the significand */
  3194.       if ((k = enormlz (yy)) > NBITS)
  3195.     ecleazs (yy);
  3196.       else
  3197.     yy[E] -= (unsigned EMUSHORT) (k - 1);
  3198.     }
  3199.   emovo (yy, y);
  3200. #endif /* not IBM */
  3201. }
  3202.  
  3203. /* Convert e-type X to IEEE 128-bit long double format E.  */
  3204.  
  3205. static void 
  3206. etoe113 (x, e)
  3207.      unsigned EMUSHORT *x, *e;
  3208. {
  3209.   unsigned EMUSHORT xi[NI];
  3210.   EMULONG exp;
  3211.   int rndsav;
  3212.  
  3213. #ifdef NANS
  3214.   if (eisnan (x))
  3215.     {
  3216.       make_nan (e, eisneg (x), TFmode);
  3217.       return;
  3218.     }
  3219. #endif
  3220.   emovi (x, xi);
  3221.   exp = (EMULONG) xi[E];
  3222. #ifdef INFINITY
  3223.   if (eisinf (x))
  3224.     goto nonorm;
  3225. #endif
  3226.   /* round off to nearest or even */
  3227.   rndsav = rndprc;
  3228.   rndprc = 113;
  3229.   emdnorm (xi, 0, 0, exp, 64);
  3230.   rndprc = rndsav;
  3231.  nonorm:
  3232.   toe113 (xi, e);
  3233. }
  3234.  
  3235. /* Convert exploded e-type X, that has already been rounded to
  3236.    113-bit precision, to IEEE 128-bit long double format Y.  */
  3237.  
  3238. static void 
  3239. toe113 (a, b)
  3240.      unsigned EMUSHORT *a, *b;
  3241. {
  3242.   register unsigned EMUSHORT *p, *q;
  3243.   unsigned EMUSHORT i;
  3244.  
  3245. #ifdef NANS
  3246.   if (eiisnan (a))
  3247.     {
  3248.       make_nan (b, eiisneg (a), TFmode);
  3249.       return;
  3250.     }
  3251. #endif
  3252.   p = a;
  3253.   if (REAL_WORDS_BIG_ENDIAN)
  3254.     q = b;
  3255.   else
  3256.     q = b + 7;            /* point to output exponent */
  3257.  
  3258.   /* If not denormal, delete the implied bit. */
  3259.   if (a[E] != 0)
  3260.     {
  3261.       eshup1 (a);
  3262.     }
  3263.   /* combine sign and exponent */
  3264.   i = *p++;
  3265.   if (REAL_WORDS_BIG_ENDIAN)
  3266.     {
  3267.       if (i)
  3268.     *q++ = *p++ | 0x8000;
  3269.       else
  3270.     *q++ = *p++;
  3271.     }
  3272.   else
  3273.     {
  3274.       if (i)
  3275.     *q-- = *p++ | 0x8000;
  3276.       else
  3277.     *q-- = *p++;
  3278.     }
  3279.   /* skip over guard word */
  3280.   ++p;
  3281.   /* move the significand */
  3282.   if (REAL_WORDS_BIG_ENDIAN)
  3283.     {
  3284.       for (i = 0; i < 7; i++)
  3285.     *q++ = *p++;
  3286.     }
  3287.   else
  3288.     {
  3289.       for (i = 0; i < 7; i++)
  3290.     *q-- = *p++;
  3291.     }
  3292. }
  3293.  
  3294. /* Convert e-type X to IEEE double extended format E.  */
  3295.  
  3296. static void 
  3297. etoe64 (x, e)
  3298.      unsigned EMUSHORT *x, *e;
  3299. {
  3300.   unsigned EMUSHORT xi[NI];
  3301.   EMULONG exp;
  3302.   int rndsav;
  3303.  
  3304. #ifdef NANS
  3305.   if (eisnan (x))
  3306.     {
  3307.       make_nan (e, eisneg (x), XFmode);
  3308.       return;
  3309.     }
  3310. #endif
  3311.   emovi (x, xi);
  3312.   /* adjust exponent for offset */
  3313.   exp = (EMULONG) xi[E];
  3314. #ifdef INFINITY
  3315.   if (eisinf (x))
  3316.     goto nonorm;
  3317. #endif
  3318.   /* round off to nearest or even */
  3319.   rndsav = rndprc;
  3320.   rndprc = 64;
  3321.   emdnorm (xi, 0, 0, exp, 64);
  3322.   rndprc = rndsav;
  3323.  nonorm:
  3324.   toe64 (xi, e);
  3325. }
  3326.  
  3327. /* Convert exploded e-type X, that has already been rounded to
  3328.    64-bit precision, to IEEE double extended format Y.  */
  3329.  
  3330. static void 
  3331. toe64 (a, b)
  3332.      unsigned EMUSHORT *a, *b;
  3333. {
  3334.   register unsigned EMUSHORT *p, *q;
  3335.   unsigned EMUSHORT i;
  3336.  
  3337. #ifdef NANS
  3338.   if (eiisnan (a))
  3339.     {
  3340.       make_nan (b, eiisneg (a), XFmode);
  3341.       return;
  3342.     }
  3343. #endif
  3344.   /* Shift denormal long double Intel format significand down one bit.  */
  3345.   if ((a[E] == 0) && ! REAL_WORDS_BIG_ENDIAN)
  3346.     eshdn1 (a);
  3347.   p = a;
  3348. #ifdef IBM
  3349.   q = b;
  3350. #endif
  3351. #ifdef DEC
  3352.   q = b + 4;
  3353. #endif
  3354. #ifdef IEEE
  3355.   if (REAL_WORDS_BIG_ENDIAN)
  3356.     q = b;
  3357.   else
  3358.     {
  3359.       q = b + 4;            /* point to output exponent */
  3360. #if LONG_DOUBLE_TYPE_SIZE == 96
  3361.       /* Clear the last two bytes of 12-byte Intel format */
  3362.       *(q+1) = 0;
  3363. #endif
  3364.     }
  3365. #endif
  3366.  
  3367.   /* combine sign and exponent */
  3368.   i = *p++;
  3369. #ifdef IBM
  3370.   if (i)
  3371.     *q++ = *p++ | 0x8000;
  3372.   else
  3373.     *q++ = *p++;
  3374.   *q++ = 0;
  3375. #endif
  3376. #ifdef DEC
  3377.   if (i)
  3378.     *q-- = *p++ | 0x8000;
  3379.   else
  3380.     *q-- = *p++;
  3381. #endif
  3382. #ifdef IEEE
  3383.   if (REAL_WORDS_BIG_ENDIAN)
  3384.     {
  3385.       if (i)
  3386.     *q++ = *p++ | 0x8000;
  3387.       else
  3388.     *q++ = *p++;
  3389.       *q++ = 0;
  3390.     }
  3391.   else
  3392.     {
  3393.       if (i)
  3394.     *q-- = *p++ | 0x8000;
  3395.       else
  3396.     *q-- = *p++;
  3397.     }
  3398. #endif
  3399.   /* skip over guard word */
  3400.   ++p;
  3401.   /* move the significand */
  3402. #ifdef IBM
  3403.   for (i = 0; i < 4; i++)
  3404.     *q++ = *p++;
  3405. #endif
  3406. #ifdef DEC
  3407.   for (i = 0; i < 4; i++)
  3408.     *q-- = *p++;
  3409. #endif
  3410. #ifdef IEEE
  3411.   if (REAL_WORDS_BIG_ENDIAN)
  3412.     {
  3413.       for (i = 0; i < 4; i++)
  3414.     *q++ = *p++;
  3415.     }
  3416.   else
  3417.     {
  3418. #ifdef INFINITY
  3419.       if (eiisinf (a))
  3420.     {
  3421.       /* Intel long double infinity significand.  */
  3422.       *q-- = 0x8000;
  3423.       *q-- = 0;
  3424.       *q-- = 0;
  3425.       *q = 0;
  3426.       return;
  3427.     }
  3428. #endif
  3429.       for (i = 0; i < 4; i++)
  3430.     *q-- = *p++;
  3431.     }
  3432. #endif
  3433. }
  3434.  
  3435. /* e type to double precision.  */
  3436.  
  3437. #ifdef DEC
  3438. /* Convert e-type X to DEC-format double E.  */
  3439.  
  3440. static void 
  3441. etoe53 (x, e)
  3442.      unsigned EMUSHORT *x, *e;
  3443. {
  3444.   etodec (x, e);        /* see etodec.c */
  3445. }
  3446.  
  3447. /* Convert exploded e-type X, that has already been rounded to
  3448.    56-bit double precision, to DEC double Y.  */
  3449.  
  3450. static void 
  3451. toe53 (x, y)
  3452.      unsigned EMUSHORT *x, *y;
  3453. {
  3454.   todec (x, y);
  3455. }
  3456.  
  3457. #else
  3458. #ifdef IBM
  3459. /* Convert e-type X to IBM 370-format double E.  */
  3460.  
  3461. static void 
  3462. etoe53 (x, e)
  3463.      unsigned EMUSHORT *x, *e;
  3464. {
  3465.   etoibm (x, e, DFmode);
  3466. }
  3467.  
  3468. /* Convert exploded e-type X, that has already been rounded to
  3469.    56-bit precision, to IBM 370 double Y.  */
  3470.  
  3471. static void 
  3472. toe53 (x, y)
  3473.      unsigned EMUSHORT *x, *y;
  3474. {
  3475.   toibm (x, y, DFmode);
  3476. }
  3477.  
  3478. #else  /* it's neither DEC nor IBM */
  3479.  
  3480. /* Convert e-type X to IEEE double E.  */
  3481.  
  3482. static void 
  3483. etoe53 (x, e)
  3484.      unsigned EMUSHORT *x, *e;
  3485. {
  3486.   unsigned EMUSHORT xi[NI];
  3487.   EMULONG exp;
  3488.   int rndsav;
  3489.  
  3490. #ifdef NANS
  3491.   if (eisnan (x))
  3492.     {
  3493.       make_nan (e, eisneg (x), DFmode);
  3494.       return;
  3495.     }
  3496. #endif
  3497.   emovi (x, xi);
  3498.   /* adjust exponent for offsets */
  3499.   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
  3500. #ifdef INFINITY
  3501.   if (eisinf (x))
  3502.     goto nonorm;
  3503. #endif
  3504.   /* round off to nearest or even */
  3505.   rndsav = rndprc;
  3506.   rndprc = 53;
  3507.   emdnorm (xi, 0, 0, exp, 64);
  3508.   rndprc = rndsav;
  3509.  nonorm:
  3510.   toe53 (xi, e);
  3511. }
  3512.  
  3513. /* Convert exploded e-type X, that has already been rounded to
  3514.    53-bit precision, to IEEE double Y.  */
  3515.  
  3516. static void 
  3517. toe53 (x, y)
  3518.      unsigned EMUSHORT *x, *y;
  3519. {
  3520.   unsigned EMUSHORT i;
  3521.   unsigned EMUSHORT *p;
  3522.  
  3523. #ifdef NANS
  3524.   if (eiisnan (x))
  3525.     {
  3526.       make_nan (y, eiisneg (x), DFmode);
  3527.       return;
  3528.     }
  3529. #endif
  3530.   p = &x[0];
  3531. #ifdef IEEE
  3532.   if (! REAL_WORDS_BIG_ENDIAN)
  3533.     y += 3;
  3534. #endif
  3535.   *y = 0;            /* output high order */
  3536.   if (*p++)
  3537.     *y = 0x8000;        /* output sign bit */
  3538.  
  3539.   i = *p++;
  3540.   if (i >= (unsigned int) 2047)
  3541.     {                /* Saturate at largest number less than infinity. */
  3542. #ifdef INFINITY
  3543.       *y |= 0x7ff0;
  3544.       if (! REAL_WORDS_BIG_ENDIAN)
  3545.     {
  3546.       *(--y) = 0;
  3547.       *(--y) = 0;
  3548.       *(--y) = 0;
  3549.     }
  3550.       else
  3551.     {
  3552.       ++y;
  3553.       *y++ = 0;
  3554.       *y++ = 0;
  3555.       *y++ = 0;
  3556.     }
  3557. #else
  3558.       *y |= (unsigned EMUSHORT) 0x7fef;
  3559.       if (! REAL_WORDS_BIG_ENDIAN)
  3560.     {
  3561.       *(--y) = 0xffff;
  3562.       *(--y) = 0xffff;
  3563.       *(--y) = 0xffff;
  3564.     }
  3565.       else
  3566.     {
  3567.       ++y;
  3568.       *y++ = 0xffff;
  3569.       *y++ = 0xffff;
  3570.       *y++ = 0xffff;
  3571.     }
  3572. #endif
  3573.       return;
  3574.     }
  3575.   if (i == 0)
  3576.     {
  3577.       eshift (x, 4);
  3578.     }
  3579.   else
  3580.     {
  3581.       i <<= 4;
  3582.       eshift (x, 5);
  3583.     }
  3584.   i |= *p++ & (unsigned EMUSHORT) 0x0f;    /* *p = xi[M] */
  3585.   *y |= (unsigned EMUSHORT) i;    /* high order output already has sign bit set */
  3586.   if (! REAL_WORDS_BIG_ENDIAN)
  3587.     {
  3588.       *(--y) = *p++;
  3589.       *(--y) = *p++;
  3590.       *(--y) = *p;
  3591.     }
  3592.   else
  3593.     {
  3594.       ++y;
  3595.       *y++ = *p++;
  3596.       *y++ = *p++;
  3597.       *y++ = *p++;
  3598.     }
  3599. }
  3600.  
  3601. #endif /* not IBM */
  3602. #endif /* not DEC */
  3603.  
  3604.  
  3605.  
  3606. /* e type to single precision.  */
  3607.  
  3608. #ifdef IBM
  3609. /* Convert e-type X to IBM 370 float E.  */
  3610.  
  3611. static void 
  3612. etoe24 (x, e)
  3613.      unsigned EMUSHORT *x, *e;
  3614. {
  3615.   etoibm (x, e, SFmode);
  3616. }
  3617.  
  3618. /* Convert exploded e-type X, that has already been rounded to
  3619.    float precision, to IBM 370 float Y.  */
  3620.  
  3621. static void 
  3622. toe24 (x, y)
  3623.      unsigned EMUSHORT *x, *y;
  3624. {
  3625.   toibm (x, y, SFmode);
  3626. }
  3627.  
  3628. #else
  3629. /* Convert e-type X to IEEE float E.  DEC float is the same as IEEE float.  */
  3630.  
  3631. static void 
  3632. etoe24 (x, e)
  3633.      unsigned EMUSHORT *x, *e;
  3634. {
  3635.   EMULONG exp;
  3636.   unsigned EMUSHORT xi[NI];
  3637.   int rndsav;
  3638.  
  3639. #ifdef NANS
  3640.   if (eisnan (x))
  3641.     {
  3642.       make_nan (e, eisneg (x), SFmode);
  3643.       return;
  3644.     }
  3645. #endif
  3646.   emovi (x, xi);
  3647.   /* adjust exponent for offsets */
  3648.   exp = (EMULONG) xi[E] - (EXONE - 0177);
  3649. #ifdef INFINITY
  3650.   if (eisinf (x))
  3651.     goto nonorm;
  3652. #endif
  3653.   /* round off to nearest or even */
  3654.   rndsav = rndprc;
  3655.   rndprc = 24;
  3656.   emdnorm (xi, 0, 0, exp, 64);
  3657.   rndprc = rndsav;
  3658.  nonorm:
  3659.   toe24 (xi, e);
  3660. }
  3661.  
  3662. /* Convert exploded e-type X, that has already been rounded to
  3663.    float precision, to IEEE float Y.  */
  3664.  
  3665. static void 
  3666. toe24 (x, y)
  3667.      unsigned EMUSHORT *x, *y;
  3668. {
  3669.   unsigned EMUSHORT i;
  3670.   unsigned EMUSHORT *p;
  3671.  
  3672. #ifdef NANS
  3673.   if (eiisnan (x))
  3674.     {
  3675.       make_nan (y, eiisneg (x), SFmode);
  3676.       return;
  3677.     }
  3678. #endif
  3679.   p = &x[0];
  3680. #ifdef IEEE
  3681.   if (! REAL_WORDS_BIG_ENDIAN)
  3682.     y += 1;
  3683. #endif
  3684. #ifdef DEC
  3685.   y += 1;
  3686. #endif
  3687.   *y = 0;            /* output high order */
  3688.   if (*p++)
  3689.     *y = 0x8000;        /* output sign bit */
  3690.  
  3691.   i = *p++;
  3692. /* Handle overflow cases. */
  3693.   if (i >= 255)
  3694.     {
  3695. #ifdef INFINITY
  3696.       *y |= (unsigned EMUSHORT) 0x7f80;
  3697. #ifdef DEC
  3698.       *(--y) = 0;
  3699. #endif
  3700. #ifdef IEEE
  3701.       if (! REAL_WORDS_BIG_ENDIAN)
  3702.     *(--y) = 0;
  3703.       else
  3704.     {
  3705.       ++y;
  3706.       *y = 0;
  3707.     }
  3708. #endif
  3709. #else  /* no INFINITY */
  3710.       *y |= (unsigned EMUSHORT) 0x7f7f;
  3711. #ifdef DEC
  3712.       *(--y) = 0xffff;
  3713. #endif
  3714. #ifdef IEEE
  3715.       if (! REAL_WORDS_BIG_ENDIAN)
  3716.     *(--y) = 0xffff;
  3717.       else
  3718.     {
  3719.       ++y;
  3720.       *y = 0xffff;
  3721.     }
  3722. #endif
  3723. #ifdef ERANGE
  3724.       errno = ERANGE;
  3725. #endif
  3726. #endif  /* no INFINITY */
  3727.       return;
  3728.     }
  3729.   if (i == 0)
  3730.     {
  3731.       eshift (x, 7);
  3732.     }
  3733.   else
  3734.     {
  3735.       i <<= 7;
  3736.       eshift (x, 8);
  3737.     }
  3738.   i |= *p++ & (unsigned EMUSHORT) 0x7f;    /* *p = xi[M] */
  3739.   /* High order output already has sign bit set.  */
  3740.   *y |= i;
  3741. #ifdef DEC
  3742.   *(--y) = *p;
  3743. #endif
  3744. #ifdef IEEE
  3745.   if (! REAL_WORDS_BIG_ENDIAN)
  3746.     *(--y) = *p;
  3747.   else
  3748.     {
  3749.       ++y;
  3750.       *y = *p;
  3751.     }
  3752. #endif
  3753. }
  3754. #endif  /* not IBM */
  3755.  
  3756. /* Compare two e type numbers. 
  3757.    Return +1 if a > b
  3758.            0 if a == b
  3759.           -1 if a < b
  3760.           -2 if either a or b is a NaN.  */
  3761.  
  3762. static int 
  3763. ecmp (a, b)
  3764.      unsigned EMUSHORT *a, *b;
  3765. {
  3766.   unsigned EMUSHORT ai[NI], bi[NI];
  3767.   register unsigned EMUSHORT *p, *q;
  3768.   register int i;
  3769.   int msign;
  3770.  
  3771. #ifdef NANS
  3772.   if (eisnan (a)  || eisnan (b))
  3773.       return (-2);
  3774. #endif
  3775.   emovi (a, ai);
  3776.   p = ai;
  3777.   emovi (b, bi);
  3778.   q = bi;
  3779.  
  3780.   if (*p != *q)
  3781.     {                /* the signs are different */
  3782.       /* -0 equals + 0 */
  3783.       for (i = 1; i < NI - 1; i++)
  3784.     {
  3785.       if (ai[i] != 0)
  3786.         goto nzro;
  3787.       if (bi[i] != 0)
  3788.         goto nzro;
  3789.     }
  3790.       return (0);
  3791.     nzro:
  3792.       if (*p == 0)
  3793.     return (1);
  3794.       else
  3795.     return (-1);
  3796.     }
  3797.   /* both are the same sign */
  3798.   if (*p == 0)
  3799.     msign = 1;
  3800.   else
  3801.     msign = -1;
  3802.   i = NI - 1;
  3803.   do
  3804.     {
  3805.       if (*p++ != *q++)
  3806.     {
  3807.       goto diff;
  3808.     }
  3809.     }
  3810.   while (--i > 0);
  3811.  
  3812.   return (0);            /* equality */
  3813.  
  3814.  diff:
  3815.  
  3816.   if (*(--p) > *(--q))
  3817.     return (msign);        /* p is bigger */
  3818.   else
  3819.     return (-msign);        /* p is littler */
  3820. }
  3821.  
  3822. /* Find e-type nearest integer to X, as floor (X + 0.5).  */
  3823.  
  3824. static void 
  3825. eround (x, y)
  3826.      unsigned EMUSHORT *x, *y;
  3827. {
  3828.   eadd (ehalf, x, y);
  3829.   efloor (y, y);
  3830. }
  3831.  
  3832. /* Convert HOST_WIDE_INT LP to e type Y.  */
  3833.  
  3834. static void 
  3835. ltoe (lp, y)
  3836.      HOST_WIDE_INT *lp;
  3837.      unsigned EMUSHORT *y;
  3838. {
  3839.   unsigned EMUSHORT yi[NI];
  3840.   unsigned HOST_WIDE_INT ll;
  3841.   int k;
  3842.  
  3843.   ecleaz (yi);
  3844.   if (*lp < 0)
  3845.     {
  3846.       /* make it positive */
  3847.       ll = (unsigned HOST_WIDE_INT) (-(*lp));
  3848.       yi[0] = 0xffff;        /* put correct sign in the e type number */
  3849.     }
  3850.   else
  3851.     {
  3852.       ll = (unsigned HOST_WIDE_INT) (*lp);
  3853.     }
  3854.   /* move the long integer to yi significand area */
  3855. #if HOST_BITS_PER_WIDE_INT == 64
  3856.   yi[M] = (unsigned EMUSHORT) (ll >> 48);
  3857.   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
  3858.   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
  3859.   yi[M + 3] = (unsigned EMUSHORT) ll;
  3860.   yi[E] = EXONE + 47;        /* exponent if normalize shift count were 0 */
  3861. #else
  3862.   yi[M] = (unsigned EMUSHORT) (ll >> 16);
  3863.   yi[M + 1] = (unsigned EMUSHORT) ll;
  3864.   yi[E] = EXONE + 15;        /* exponent if normalize shift count were 0 */
  3865. #endif
  3866.  
  3867.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  3868.     ecleaz (yi);        /* it was zero */
  3869.   else
  3870.     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
  3871.   emovo (yi, y);        /* output the answer */
  3872. }
  3873.  
  3874. /* Convert unsigned HOST_WIDE_INT LP to e type Y.  */
  3875.  
  3876. static void 
  3877. ultoe (lp, y)
  3878.      unsigned HOST_WIDE_INT *lp;
  3879.      unsigned EMUSHORT *y;
  3880. {
  3881.   unsigned EMUSHORT yi[NI];
  3882.   unsigned HOST_WIDE_INT ll;
  3883.   int k;
  3884.  
  3885.   ecleaz (yi);
  3886.   ll = *lp;
  3887.  
  3888.   /* move the long integer to ayi significand area */
  3889. #if HOST_BITS_PER_WIDE_INT == 64
  3890.   yi[M] = (unsigned EMUSHORT) (ll >> 48);
  3891.   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
  3892.   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
  3893.   yi[M + 3] = (unsigned EMUSHORT) ll;
  3894.   yi[E] = EXONE + 47;        /* exponent if normalize shift count were 0 */
  3895. #else
  3896.   yi[M] = (unsigned EMUSHORT) (ll >> 16);
  3897.   yi[M + 1] = (unsigned EMUSHORT) ll;
  3898.   yi[E] = EXONE + 15;        /* exponent if normalize shift count were 0 */
  3899. #endif
  3900.  
  3901.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  3902.     ecleaz (yi);        /* it was zero */
  3903.   else
  3904.     yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
  3905.   emovo (yi, y);        /* output the answer */
  3906. }
  3907.  
  3908.  
  3909. /* Find signed HOST_WIDE_INT integer I and floating point fractional
  3910.    part FRAC of e-type (packed internal format) floating point input X.
  3911.    The integer output I has the sign of the input, except that
  3912.    positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
  3913.    The output e-type fraction FRAC is the positive fractional
  3914.    part of abs (X).  */
  3915.  
  3916. static void 
  3917. eifrac (x, i, frac)
  3918.      unsigned EMUSHORT *x;
  3919.      HOST_WIDE_INT *i;
  3920.      unsigned EMUSHORT *frac;
  3921. {
  3922.   unsigned EMUSHORT xi[NI];
  3923.   int j, k;
  3924.   unsigned HOST_WIDE_INT ll;
  3925.  
  3926.   emovi (x, xi);
  3927.   k = (int) xi[E] - (EXONE - 1);
  3928.   if (k <= 0)
  3929.     {
  3930.       /* if exponent <= 0, integer = 0 and real output is fraction */
  3931.       *i = 0L;
  3932.       emovo (xi, frac);
  3933.       return;
  3934.     }
  3935.   if (k > (HOST_BITS_PER_WIDE_INT - 1))
  3936.     {
  3937.       /* long integer overflow: output large integer
  3938.      and correct fraction  */
  3939.       if (xi[0])
  3940.     *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
  3941.       else
  3942.     {
  3943. #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
  3944.       /* In this case, let it overflow and convert as if unsigned.  */
  3945.       euifrac (x, &ll, frac);
  3946.       *i = (HOST_WIDE_INT) ll;
  3947.       return;
  3948. #else
  3949.       /* In other cases, return the largest positive integer.  */
  3950.       *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
  3951. #endif
  3952.     }
  3953.       eshift (xi, k);
  3954.       if (extra_warnings)
  3955.     warning ("overflow on truncation to integer");
  3956.     }
  3957.   else if (k > 16)
  3958.     {
  3959.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  3960.      then shift up by 16's.  */
  3961.       j = k - ((k >> 4) << 4);
  3962.       eshift (xi, j);
  3963.       ll = xi[M];
  3964.       k -= j;
  3965.       do
  3966.     {
  3967.       eshup6 (xi);
  3968.       ll = (ll << 16) | xi[M];
  3969.     }
  3970.       while ((k -= 16) > 0);
  3971.       *i = ll;
  3972.       if (xi[0])
  3973.     *i = -(*i);
  3974.     }
  3975.   else
  3976.       {
  3977.         /* shift not more than 16 bits */
  3978.           eshift (xi, k);
  3979.         *i = (HOST_WIDE_INT) xi[M] & 0xffff;
  3980.         if (xi[0])
  3981.       *i = -(*i);
  3982.       }
  3983.   xi[0] = 0;
  3984.   xi[E] = EXONE - 1;
  3985.   xi[M] = 0;
  3986.   if ((k = enormlz (xi)) > NBITS)
  3987.     ecleaz (xi);
  3988.   else
  3989.     xi[E] -= (unsigned EMUSHORT) k;
  3990.  
  3991.   emovo (xi, frac);
  3992. }
  3993.  
  3994.  
  3995. /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
  3996.    FRAC of e-type X.  A negative input yields integer output = 0 but
  3997.    correct fraction.  */
  3998.  
  3999. static void 
  4000. euifrac (x, i, frac)
  4001.      unsigned EMUSHORT *x;
  4002.      unsigned HOST_WIDE_INT *i;
  4003.      unsigned EMUSHORT *frac;
  4004. {
  4005.   unsigned HOST_WIDE_INT ll;
  4006.   unsigned EMUSHORT xi[NI];
  4007.   int j, k;
  4008.  
  4009.   emovi (x, xi);
  4010.   k = (int) xi[E] - (EXONE - 1);
  4011.   if (k <= 0)
  4012.     {
  4013.       /* if exponent <= 0, integer = 0 and argument is fraction */
  4014.       *i = 0L;
  4015.       emovo (xi, frac);
  4016.       return;
  4017.     }
  4018.   if (k > HOST_BITS_PER_WIDE_INT)
  4019.     {
  4020.       /* Long integer overflow: output large integer
  4021.      and correct fraction.
  4022.      Note, the BSD microvax compiler says that ~(0UL)
  4023.      is a syntax error.  */
  4024.       *i = ~(0L);
  4025.       eshift (xi, k);
  4026.       if (extra_warnings)
  4027.     warning ("overflow on truncation to unsigned integer");
  4028.     }
  4029.   else if (k > 16)
  4030.     {
  4031.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  4032.      then shift up by 16's.  */
  4033.       j = k - ((k >> 4) << 4);
  4034.       eshift (xi, j);
  4035.       ll = xi[M];
  4036.       k -= j;
  4037.       do
  4038.     {
  4039.       eshup6 (xi);
  4040.       ll = (ll << 16) | xi[M];
  4041.     }
  4042.       while ((k -= 16) > 0);
  4043.       *i = ll;
  4044.     }
  4045.   else
  4046.     {
  4047.       /* shift not more than 16 bits */
  4048.       eshift (xi, k);
  4049.       *i = (HOST_WIDE_INT) xi[M] & 0xffff;
  4050.     }
  4051.  
  4052.   if (xi[0])  /* A negative value yields unsigned integer 0. */
  4053.     *i = 0L;
  4054.  
  4055.   xi[0] = 0;
  4056.   xi[E] = EXONE - 1;
  4057.   xi[M] = 0;
  4058.   if ((k = enormlz (xi)) > NBITS)
  4059.     ecleaz (xi);
  4060.   else
  4061.     xi[E] -= (unsigned EMUSHORT) k;
  4062.  
  4063.   emovo (xi, frac);
  4064. }
  4065.  
  4066. /* Shift the significand of exploded e-type X up or down by SC bits.  */
  4067.  
  4068. static int 
  4069. eshift (x, sc)
  4070.      unsigned EMUSHORT *x;
  4071.      int sc;
  4072. {
  4073.   unsigned EMUSHORT lost;
  4074.   unsigned EMUSHORT *p;
  4075.  
  4076.   if (sc == 0)
  4077.     return (0);
  4078.  
  4079.   lost = 0;
  4080.   p = x + NI - 1;
  4081.  
  4082.   if (sc < 0)
  4083.     {
  4084.       sc = -sc;
  4085.       while (sc >= 16)
  4086.     {
  4087.       lost |= *p;        /* remember lost bits */
  4088.       eshdn6 (x);
  4089.       sc -= 16;
  4090.     }
  4091.  
  4092.       while (sc >= 8)
  4093.     {
  4094.       lost |= *p & 0xff;
  4095.       eshdn8 (x);
  4096.       sc -= 8;
  4097.     }
  4098.  
  4099.       while (sc > 0)
  4100.     {
  4101.       lost |= *p & 1;
  4102.       eshdn1 (x);
  4103.       sc -= 1;
  4104.     }
  4105.     }
  4106.   else
  4107.     {
  4108.       while (sc >= 16)
  4109.     {
  4110.       eshup6 (x);
  4111.       sc -= 16;
  4112.     }
  4113.  
  4114.       while (sc >= 8)
  4115.     {
  4116.       eshup8 (x);
  4117.       sc -= 8;
  4118.     }
  4119.  
  4120.       while (sc > 0)
  4121.     {
  4122.       eshup1 (x);
  4123.       sc -= 1;
  4124.     }
  4125.     }
  4126.   if (lost)
  4127.     lost = 1;
  4128.   return ((int) lost);
  4129. }
  4130.  
  4131. /* Shift normalize the significand area of exploded e-type X.
  4132.    Return the shift count (up = positive).  */
  4133.  
  4134. static int 
  4135. enormlz (x)
  4136.      unsigned EMUSHORT x[];
  4137. {
  4138.   register unsigned EMUSHORT *p;
  4139.   int sc;
  4140.  
  4141.   sc = 0;
  4142.   p = &x[M];
  4143.   if (*p != 0)
  4144.     goto normdn;
  4145.   ++p;
  4146.   if (*p & 0x8000)
  4147.     return (0);            /* already normalized */
  4148.   while (*p == 0)
  4149.     {
  4150.       eshup6 (x);
  4151.       sc += 16;
  4152.  
  4153.       /* With guard word, there are NBITS+16 bits available.
  4154.        Return true if all are zero.  */
  4155.       if (sc > NBITS)
  4156.     return (sc);
  4157.     }
  4158.   /* see if high byte is zero */
  4159.   while ((*p & 0xff00) == 0)
  4160.     {
  4161.       eshup8 (x);
  4162.       sc += 8;
  4163.     }
  4164.   /* now shift 1 bit at a time */
  4165.   while ((*p & 0x8000) == 0)
  4166.     {
  4167.       eshup1 (x);
  4168.       sc += 1;
  4169.       if (sc > NBITS)
  4170.     {
  4171.       mtherr ("enormlz", UNDERFLOW);
  4172.       return (sc);
  4173.     }
  4174.     }
  4175.   return (sc);
  4176.  
  4177.   /* Normalize by shifting down out of the high guard word
  4178.      of the significand */
  4179.  normdn:
  4180.  
  4181.   if (*p & 0xff00)
  4182.     {
  4183.       eshdn8 (x);
  4184.       sc -= 8;
  4185.     }
  4186.   while (*p != 0)
  4187.     {
  4188.       eshdn1 (x);
  4189.       sc -= 1;
  4190.  
  4191.       if (sc < -NBITS)
  4192.     {
  4193.       mtherr ("enormlz", OVERFLOW);
  4194.       return (sc);
  4195.     }
  4196.     }
  4197.   return (sc);
  4198. }
  4199.  
  4200. /* Powers of ten used in decimal <-> binary conversions.  */
  4201.  
  4202. #define NTEN 12
  4203. #define MAXP 4096
  4204.  
  4205. #if LONG_DOUBLE_TYPE_SIZE == 128
  4206. static unsigned EMUSHORT etens[NTEN + 1][NE] =
  4207. {
  4208.   {0x6576, 0x4a92, 0x804a, 0x153f,
  4209.    0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
  4210.   {0x6a32, 0xce52, 0x329a, 0x28ce,
  4211.    0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
  4212.   {0x526c, 0x50ce, 0xf18b, 0x3d28,
  4213.    0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
  4214.   {0x9c66, 0x58f8, 0xbc50, 0x5c54,
  4215.    0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
  4216.   {0x851e, 0xeab7, 0x98fe, 0x901b,
  4217.    0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
  4218.   {0x0235, 0x0137, 0x36b1, 0x336c,
  4219.    0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
  4220.   {0x50f8, 0x25fb, 0xc76b, 0x6b71,
  4221.    0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
  4222.   {0x0000, 0x0000, 0x0000, 0x0000,
  4223.    0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
  4224.   {0x0000, 0x0000, 0x0000, 0x0000,
  4225.    0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
  4226.   {0x0000, 0x0000, 0x0000, 0x0000,
  4227.    0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
  4228.   {0x0000, 0x0000, 0x0000, 0x0000,
  4229.    0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
  4230.   {0x0000, 0x0000, 0x0000, 0x0000,
  4231.    0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
  4232.   {0x0000, 0x0000, 0x0000, 0x0000,
  4233.    0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
  4234. };
  4235.  
  4236. static unsigned EMUSHORT emtens[NTEN + 1][NE] =
  4237. {
  4238.   {0x2030, 0xcffc, 0xa1c3, 0x8123,
  4239.    0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
  4240.   {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
  4241.    0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
  4242.   {0xf53f, 0xf698, 0x6bd3, 0x0158,
  4243.    0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
  4244.   {0xe731, 0x04d4, 0xe3f2, 0xd332,
  4245.    0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
  4246.   {0xa23e, 0x5308, 0xfefb, 0x1155,
  4247.    0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
  4248.   {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
  4249.    0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
  4250.   {0x2a20, 0x6224, 0x47b3, 0x98d7,
  4251.    0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
  4252.   {0x0b5b, 0x4af2, 0xa581, 0x18ed,
  4253.    0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
  4254.   {0xbf71, 0xa9b3, 0x7989, 0xbe68,
  4255.    0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
  4256.   {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
  4257.    0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
  4258.   {0xc155, 0xa4a8, 0x404e, 0x6113,
  4259.    0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
  4260.   {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
  4261.    0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
  4262.   {0xcccd, 0xcccc, 0xcccc, 0xcccc,
  4263.    0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
  4264. };
  4265. #else
  4266. /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
  4267. static unsigned EMUSHORT etens[NTEN + 1][NE] =
  4268. {
  4269.   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
  4270.   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
  4271.   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
  4272.   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
  4273.   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
  4274.   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
  4275.   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
  4276.   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
  4277.   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
  4278.   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
  4279.   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
  4280.   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
  4281.   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
  4282. };
  4283.  
  4284. static unsigned EMUSHORT emtens[NTEN + 1][NE] =
  4285. {
  4286.   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
  4287.   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
  4288.   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
  4289.   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
  4290.   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
  4291.   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
  4292.   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
  4293.   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
  4294.   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
  4295.   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
  4296.   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
  4297.   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
  4298.   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
  4299. };
  4300. #endif
  4301.  
  4302. /* Convert float value X to ASCII string STRING with NDIG digits after
  4303.    the decimal point.  */
  4304.  
  4305. static void 
  4306. e24toasc (x, string, ndigs)
  4307.      unsigned EMUSHORT x[];
  4308.      char *string;
  4309.      int ndigs;
  4310. {
  4311.   unsigned EMUSHORT w[NI];
  4312.  
  4313.   e24toe (x, w);
  4314.   etoasc (w, string, ndigs);
  4315. }
  4316.  
  4317. /* Convert double value X to ASCII string STRING with NDIG digits after
  4318.    the decimal point.  */
  4319.  
  4320. static void 
  4321. e53toasc (x, string, ndigs)
  4322.      unsigned EMUSHORT x[];
  4323.      char *string;
  4324.      int ndigs;
  4325. {
  4326.   unsigned EMUSHORT w[NI];
  4327.  
  4328.   e53toe (x, w);
  4329.   etoasc (w, string, ndigs);
  4330. }
  4331.  
  4332. /* Convert double extended value X to ASCII string STRING with NDIG digits
  4333.    after the decimal point.  */
  4334.  
  4335. static void 
  4336. e64toasc (x, string, ndigs)
  4337.      unsigned EMUSHORT x[];
  4338.      char *string;
  4339.      int ndigs;
  4340. {
  4341.   unsigned EMUSHORT w[NI];
  4342.  
  4343.   e64toe (x, w);
  4344.   etoasc (w, string, ndigs);
  4345. }
  4346.  
  4347. /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
  4348.    after the decimal point.  */
  4349.  
  4350. static void 
  4351. e113toasc (x, string, ndigs)
  4352.      unsigned EMUSHORT x[];
  4353.      char *string;
  4354.      int ndigs;
  4355. {
  4356.   unsigned EMUSHORT w[NI];
  4357.  
  4358.   e113toe (x, w);
  4359.   etoasc (w, string, ndigs);
  4360. }
  4361.  
  4362. /* Convert e-type X to ASCII string STRING with NDIGS digits after
  4363.    the decimal point.  */
  4364.  
  4365. static char wstring[80];    /* working storage for ASCII output */
  4366.  
  4367. static void 
  4368. etoasc (x, string, ndigs)
  4369.      unsigned EMUSHORT x[];
  4370.      char *string;
  4371.      int ndigs;
  4372. {
  4373.   EMUSHORT digit;
  4374.   unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
  4375.   unsigned EMUSHORT *p, *r, *ten;
  4376.   unsigned EMUSHORT sign;
  4377.   int i, j, k, expon, rndsav;
  4378.   char *s, *ss;
  4379.   unsigned EMUSHORT m;
  4380.  
  4381.  
  4382.   rndsav = rndprc;
  4383.   ss = string;
  4384.   s = wstring;
  4385.   *ss = '\0';
  4386.   *s = '\0';
  4387. #ifdef NANS
  4388.   if (eisnan (x))
  4389.     {
  4390.       sprintf (wstring, " NaN ");
  4391.       goto bxit;
  4392.     }
  4393. #endif
  4394.   rndprc = NBITS;        /* set to full precision */
  4395.   emov (x, y);            /* retain external format */
  4396.   if (y[NE - 1] & 0x8000)
  4397.     {
  4398.       sign = 0xffff;
  4399.       y[NE - 1] &= 0x7fff;
  4400.     }
  4401.   else
  4402.     {
  4403.       sign = 0;
  4404.     }
  4405.   expon = 0;
  4406.   ten = &etens[NTEN][0];
  4407.   emov (eone, t);
  4408.   /* Test for zero exponent */
  4409.   if (y[NE - 1] == 0)
  4410.     {
  4411.       for (k = 0; k < NE - 1; k++)
  4412.     {
  4413.       if (y[k] != 0)
  4414.         goto tnzro;        /* denormalized number */
  4415.     }
  4416.       goto isone;        /* valid all zeros */
  4417.     }
  4418.  tnzro:
  4419.  
  4420.   /* Test for infinity. */
  4421.   if (y[NE - 1] == 0x7fff)
  4422.     {
  4423.       if (sign)
  4424. #ifdef __amigados__
  4425.     sprintf (wstring, " -Infinity ");
  4426.       else
  4427.     sprintf (wstring, " Infinity ");
  4428. #else
  4429.     sprintf (wstring, " -NaN ");
  4430.       else
  4431.     sprintf (wstring, " NaN ");
  4432. #endif
  4433.       goto bxit;
  4434.     }
  4435.  
  4436.   /* Test for exponent nonzero but significand denormalized.
  4437.    * This is an error condition.
  4438.    */
  4439.   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
  4440.     {
  4441.       mtherr ("etoasc", DOMAIN);
  4442.       sprintf (wstring, "NaN");
  4443.       goto bxit;
  4444.     }
  4445.  
  4446.   /* Compare to 1.0 */
  4447.   i = ecmp (eone, y);
  4448.   if (i == 0)
  4449.     goto isone;
  4450.  
  4451.   if (i == -2)
  4452.     abort ();
  4453.  
  4454.   if (i < 0)
  4455.     {                /* Number is greater than 1 */
  4456.       /* Convert significand to an integer and strip trailing decimal zeros. */
  4457.       emov (y, u);
  4458.       u[NE - 1] = EXONE + NBITS - 1;
  4459.  
  4460.       p = &etens[NTEN - 4][0];
  4461.       m = 16;
  4462.       do
  4463.     {
  4464.       ediv (p, u, t);
  4465.       efloor (t, w);
  4466.       for (j = 0; j < NE - 1; j++)
  4467.         {
  4468.           if (t[j] != w[j])
  4469.         goto noint;
  4470.         }
  4471.       emov (t, u);
  4472.       expon += (int) m;
  4473.     noint:
  4474.       p += NE;
  4475.       m >>= 1;
  4476.     }
  4477.       while (m != 0);
  4478.  
  4479.       /* Rescale from integer significand */
  4480.       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
  4481.       emov (u, y);
  4482.       /* Find power of 10 */
  4483.       emov (eone, t);
  4484.       m = MAXP;
  4485.       p = &etens[0][0];
  4486.       /* An unordered compare result shouldn't happen here. */
  4487.       while (ecmp (ten, u) <= 0)
  4488.     {
  4489.       if (ecmp (p, u) <= 0)
  4490.         {
  4491.           ediv (p, u, u);
  4492.           emul (p, t, t);
  4493.           expon += (int) m;
  4494.         }
  4495.       m >>= 1;
  4496.       if (m == 0)
  4497.         break;
  4498.       p += NE;
  4499.     }
  4500.     }
  4501.   else
  4502.     {                /* Number is less than 1.0 */
  4503.       /* Pad significand with trailing decimal zeros. */
  4504.       if (y[NE - 1] == 0)
  4505.     {
  4506.       while ((y[NE - 2] & 0x8000) == 0)
  4507.         {
  4508.           emul (ten, y, y);
  4509.           expon -= 1;
  4510.         }
  4511.     }
  4512.       else
  4513.     {
  4514.       emovi (y, w);
  4515.       for (i = 0; i < NDEC + 1; i++)
  4516.         {
  4517.           if ((w[NI - 1] & 0x7) != 0)
  4518.         break;
  4519.           /* multiply by 10 */
  4520.           emovz (w, u);
  4521.           eshdn1 (u);
  4522.           eshdn1 (u);
  4523.           eaddm (w, u);
  4524.           u[1] += 3;
  4525.           while (u[2] != 0)
  4526.         {
  4527.           eshdn1 (u);
  4528.           u[1] += 1;
  4529.         }
  4530.           if (u[NI - 1] != 0)
  4531.         break;
  4532.           if (eone[NE - 1] <= u[1])
  4533.         break;
  4534.           emovz (u, w);
  4535.           expon -= 1;
  4536.         }
  4537.       emovo (w, y);
  4538.     }
  4539.       k = -MAXP;
  4540.       p = &emtens[0][0];
  4541.       r = &etens[0][0];
  4542.       emov (y, w);
  4543.       emov (eone, t);
  4544.       while (ecmp (eone, w) > 0)
  4545.     {
  4546.       if (ecmp (p, w) >= 0)
  4547.         {
  4548.           emul (r, w, w);
  4549.           emul (r, t, t);
  4550.           expon += k;
  4551.         }
  4552.       k /= 2;
  4553.       if (k == 0)
  4554.         break;
  4555.       p += NE;
  4556.       r += NE;
  4557.     }
  4558.       ediv (t, eone, t);
  4559.     }
  4560.  isone:
  4561.   /* Find the first (leading) digit. */
  4562.   emovi (t, w);
  4563.   emovz (w, t);
  4564.   emovi (y, w);
  4565.   emovz (w, y);
  4566.   eiremain (t, y);
  4567.   digit = equot[NI - 1];
  4568.   while ((digit == 0) && (ecmp (y, ezero) != 0))
  4569.     {
  4570.       eshup1 (y);
  4571.       emovz (y, u);
  4572.       eshup1 (u);
  4573.       eshup1 (u);
  4574.       eaddm (u, y);
  4575.       eiremain (t, y);
  4576.       digit = equot[NI - 1];
  4577.       expon -= 1;
  4578.     }
  4579.   s = wstring;
  4580.   if (sign)
  4581.     *s++ = '-';
  4582.   else
  4583.     *s++ = ' ';
  4584.   /* Examine number of digits requested by caller. */
  4585.   if (ndigs < 0)
  4586.     ndigs = 0;
  4587.   if (ndigs > NDEC)
  4588.     ndigs = NDEC;
  4589.   if (digit == 10)
  4590.     {
  4591.       *s++ = '1';
  4592.       *s++ = '.';
  4593.       if (ndigs > 0)
  4594.     {
  4595.       *s++ = '0';
  4596.       ndigs -= 1;
  4597.     }
  4598.       expon += 1;
  4599.     }
  4600.   else
  4601.     {
  4602.       *s++ = (char)digit + '0';
  4603.       *s++ = '.';
  4604.     }
  4605.   /* Generate digits after the decimal point. */
  4606.   for (k = 0; k <= ndigs; k++)
  4607.     {
  4608.       /* multiply current number by 10, without normalizing */
  4609.       eshup1 (y);
  4610.       emovz (y, u);
  4611.       eshup1 (u);
  4612.       eshup1 (u);
  4613.       eaddm (u, y);
  4614.       eiremain (t, y);
  4615.       *s++ = (char) equot[NI - 1] + '0';
  4616.     }
  4617.   digit = equot[NI - 1];
  4618.   --s;
  4619.   ss = s;
  4620.   /* round off the ASCII string */
  4621.   if (digit > 4)
  4622.     {
  4623.       /* Test for critical rounding case in ASCII output. */
  4624.       if (digit == 5)
  4625.     {
  4626.       emovo (y, t);
  4627.       if (ecmp (t, ezero) != 0)
  4628.         goto roun;        /* round to nearest */
  4629.       if ((*(s - 1) & 1) == 0)
  4630.         goto doexp;        /* round to even */
  4631.     }
  4632.       /* Round up and propagate carry-outs */
  4633.     roun:
  4634.       --s;
  4635.       k = *s & 0x7f;
  4636.       /* Carry out to most significant digit? */
  4637.       if (k == '.')
  4638.     {
  4639.       --s;
  4640.       k = *s;
  4641.       k += 1;
  4642.       *s = (char) k;
  4643.       /* Most significant digit carries to 10? */
  4644.       if (k > '9')
  4645.         {
  4646.           expon += 1;
  4647.           *s = '1';
  4648.         }
  4649.       goto doexp;
  4650.     }
  4651.       /* Round up and carry out from less significant digits */
  4652.       k += 1;
  4653.       *s = (char) k;
  4654.       if (k > '9')
  4655.     {
  4656.       *s = '0';
  4657.       goto roun;
  4658.     }
  4659.     }
  4660.  doexp:
  4661.   /*
  4662.      if (expon >= 0)
  4663.      sprintf (ss, "e+%d", expon);
  4664.      else
  4665.      sprintf (ss, "e%d", expon);
  4666.      */
  4667.   sprintf (ss, "e%d", expon);
  4668.  bxit:
  4669.   rndprc = rndsav;
  4670.   /* copy out the working string */
  4671.   s = string;
  4672.   ss = wstring;
  4673.   while (*ss == ' ')        /* strip possible leading space */
  4674.     ++ss;
  4675.   while ((*s++ = *ss++) != '\0')
  4676.     ;
  4677. }
  4678.  
  4679.  
  4680. /* Convert ASCII string to floating point.
  4681.  
  4682.    Numeric input is a free format decimal number of any length, with
  4683.    or without decimal point.  Entering E after the number followed by an
  4684.    integer number causes the second number to be interpreted as a power of
  4685.    10 to be multiplied by the first number (i.e., "scientific" notation).  */
  4686.  
  4687. /* Convert ASCII string S to single precision float value Y.  */
  4688.  
  4689. static void 
  4690. asctoe24 (s, y)
  4691.      char *s;
  4692.      unsigned EMUSHORT *y;
  4693. {
  4694.   asctoeg (s, y, 24);
  4695. }
  4696.  
  4697.  
  4698. /* Convert ASCII string S to double precision value Y.  */
  4699.  
  4700. static void 
  4701. asctoe53 (s, y)
  4702.      char *s;
  4703.      unsigned EMUSHORT *y;
  4704. {
  4705. #if defined(DEC) || defined(IBM)
  4706.   asctoeg (s, y, 56);
  4707. #else
  4708.   asctoeg (s, y, 53);
  4709. #endif
  4710. }
  4711.  
  4712.  
  4713. /* Convert ASCII string S to double extended value Y.  */
  4714.  
  4715. static void 
  4716. asctoe64 (s, y)
  4717.      char *s;
  4718.      unsigned EMUSHORT *y;
  4719. {
  4720.   asctoeg (s, y, 64);
  4721. }
  4722.  
  4723. /* Convert ASCII string S to 128-bit long double Y.  */
  4724.  
  4725. static void 
  4726. asctoe113 (s, y)
  4727.      char *s;
  4728.      unsigned EMUSHORT *y;
  4729. {
  4730.   asctoeg (s, y, 113);
  4731. }
  4732.  
  4733. /* Convert ASCII string S to e type Y.  */
  4734.  
  4735. static void 
  4736. asctoe (s, y)
  4737.      char *s;
  4738.      unsigned EMUSHORT *y;
  4739. {
  4740.   asctoeg (s, y, NBITS);
  4741. }
  4742.  
  4743. /* Convert ASCII string SS to e type Y, with a specified rounding precision
  4744.    of OPREC bits. */
  4745.  
  4746. static void 
  4747. asctoeg (ss, y, oprec)
  4748.      char *ss;
  4749.      unsigned EMUSHORT *y;
  4750.      int oprec;
  4751. {
  4752.   unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
  4753.   int esign, decflg, sgnflg, nexp, exp, prec, lost;
  4754.   int k, trail, c, rndsav;
  4755.   EMULONG lexp;
  4756.   unsigned EMUSHORT nsign, *p;
  4757.   char *sp, *s, *lstr;
  4758.  
  4759.   /* Copy the input string. */
  4760.   lstr = (char *) alloca (strlen (ss) + 1);
  4761.   s = ss;
  4762.   while (*s == ' ')        /* skip leading spaces */
  4763.     ++s;
  4764.   sp = lstr;
  4765.   while ((*sp++ = *s++) != '\0')
  4766.     ;
  4767.   s = lstr;
  4768.  
  4769.   rndsav = rndprc;
  4770.   rndprc = NBITS;        /* Set to full precision */
  4771.   lost = 0;
  4772.   nsign = 0;
  4773.   decflg = 0;
  4774.   sgnflg = 0;
  4775.   nexp = 0;
  4776.   exp = 0;
  4777.   prec = 0;
  4778.   ecleaz (yy);
  4779.   trail = 0;
  4780.  
  4781.  nxtcom:
  4782.   k = *s - '0';
  4783.   if ((k >= 0) && (k <= 9))
  4784.     {
  4785.       /* Ignore leading zeros */
  4786.       if ((prec == 0) && (decflg == 0) && (k == 0))
  4787.     goto donchr;
  4788.       /* Identify and strip trailing zeros after the decimal point. */
  4789.       if ((trail == 0) && (decflg != 0))
  4790.     {
  4791.       sp = s;
  4792.       while ((*sp >= '0') && (*sp <= '9'))
  4793.         ++sp;
  4794.       /* Check for syntax error */
  4795.       c = *sp & 0x7f;
  4796.       if ((c != 'e') && (c != 'E') && (c != '\0')
  4797.           && (c != '\n') && (c != '\r') && (c != ' ')
  4798.           && (c != ','))
  4799.         goto error;
  4800.       --sp;
  4801.       while (*sp == '0')
  4802.         *sp-- = 'z';
  4803.       trail = 1;
  4804.       if (*s == 'z')
  4805.         goto donchr;
  4806.     }
  4807.  
  4808.       /* If enough digits were given to more than fill up the yy register,
  4809.      continuing until overflow into the high guard word yy[2]
  4810.      guarantees that there will be a roundoff bit at the top
  4811.      of the low guard word after normalization.  */
  4812.  
  4813.       if (yy[2] == 0)
  4814.     {
  4815.       if (decflg)
  4816.         nexp += 1;        /* count digits after decimal point */
  4817.       eshup1 (yy);        /* multiply current number by 10 */
  4818.       emovz (yy, xt);
  4819.       eshup1 (xt);
  4820.       eshup1 (xt);
  4821.       eaddm (xt, yy);
  4822.       ecleaz (xt);
  4823.       xt[NI - 2] = (unsigned EMUSHORT) k;
  4824.       eaddm (xt, yy);
  4825.     }
  4826.       else
  4827.     {
  4828.       /* Mark any lost non-zero digit.  */
  4829.       lost |= k;
  4830.       /* Count lost digits before the decimal point.  */
  4831.       if (decflg == 0)
  4832.         nexp -= 1;
  4833.     }
  4834.       prec += 1;
  4835.       goto donchr;
  4836.     }
  4837.  
  4838.   switch (*s)
  4839.     {
  4840.     case 'z':
  4841.       break;
  4842.     case 'E':
  4843.     case 'e':
  4844.       goto expnt;
  4845.     case '.':            /* decimal point */
  4846.       if (decflg)
  4847.     goto error;
  4848.       ++decflg;
  4849.       break;
  4850.     case '-':
  4851.       nsign = 0xffff;
  4852.       if (sgnflg)
  4853.     goto error;
  4854.       ++sgnflg;
  4855.       break;
  4856.     case '+':
  4857.       if (sgnflg)
  4858.     goto error;
  4859.       ++sgnflg;
  4860.       break;
  4861.     case ',':
  4862.     case ' ':
  4863.     case '\0':
  4864.     case '\n':
  4865.     case '\r':
  4866.       goto daldone;
  4867.     case 'i':
  4868.     case 'I':
  4869.       goto infinite;
  4870.     default:
  4871.     error:
  4872. #ifdef NANS
  4873.       einan (yy);
  4874. #else
  4875.       mtherr ("asctoe", DOMAIN);
  4876.       eclear (yy);
  4877. #endif
  4878.       goto aexit;
  4879.     }
  4880.  donchr:
  4881.   ++s;
  4882.   goto nxtcom;
  4883.  
  4884.   /* Exponent interpretation */
  4885.  expnt:
  4886.  
  4887.   esign = 1;
  4888.   exp = 0;
  4889.   ++s;
  4890.   /* check for + or - */
  4891.   if (*s == '-')
  4892.     {
  4893.       esign = -1;
  4894.       ++s;
  4895.     }
  4896.   if (*s == '+')
  4897.     ++s;
  4898.   while ((*s >= '0') && (*s <= '9'))
  4899.     {
  4900.       exp *= 10;
  4901.       exp += *s++ - '0';
  4902.       if (exp > -(MINDECEXP))
  4903.     {
  4904.       if (esign < 0)
  4905.         goto zero;
  4906.       else
  4907.         goto infinite;
  4908.     }
  4909.     }
  4910.   if (esign < 0)
  4911.     exp = -exp;
  4912.   if (exp > MAXDECEXP)
  4913.     {
  4914.  infinite:
  4915.       ecleaz (yy);
  4916.       yy[E] = 0x7fff;        /* infinity */
  4917.       goto aexit;
  4918.     }
  4919.   if (exp < MINDECEXP)
  4920.     {
  4921.  zero:
  4922.       ecleaz (yy);
  4923.       goto aexit;
  4924.     }
  4925.  
  4926.  daldone:
  4927.   nexp = exp - nexp;
  4928.   /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
  4929.   while ((nexp > 0) && (yy[2] == 0))
  4930.     {
  4931.       emovz (yy, xt);
  4932.       eshup1 (xt);
  4933.       eshup1 (xt);
  4934.       eaddm (yy, xt);
  4935.       eshup1 (xt);
  4936.       if (xt[2] != 0)
  4937.     break;
  4938.       nexp -= 1;
  4939.       emovz (xt, yy);
  4940.     }
  4941.   if ((k = enormlz (yy)) > NBITS)
  4942.     {
  4943.       ecleaz (yy);
  4944.       goto aexit;
  4945.     }
  4946.   lexp = (EXONE - 1 + NBITS) - k;
  4947.   emdnorm (yy, lost, 0, lexp, 64);
  4948.  
  4949.   /* Convert to external format:
  4950.  
  4951.      Multiply by 10**nexp.  If precision is 64 bits,
  4952.      the maximum relative error incurred in forming 10**n
  4953.      for 0 <= n <= 324 is 8.2e-20, at 10**180.
  4954.      For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
  4955.      For 0 >= n >= -999, it is -1.55e-19 at 10**-435.  */
  4956.  
  4957.   lexp = yy[E];
  4958.   if (nexp == 0)
  4959.     {
  4960.       k = 0;
  4961.       goto expdon;
  4962.     }
  4963.   esign = 1;
  4964.   if (nexp < 0)
  4965.     {
  4966.       nexp = -nexp;
  4967.       esign = -1;
  4968.       if (nexp > 4096)
  4969.     {
  4970.       /* Punt.  Can't handle this without 2 divides. */
  4971.       emovi (etens[0], tt);
  4972.       lexp -= tt[E];
  4973.       k = edivm (tt, yy);
  4974.       lexp += EXONE;
  4975.       nexp -= 4096;
  4976.     }
  4977.     }
  4978.   p = &etens[NTEN][0];
  4979.   emov (eone, xt);
  4980.   exp = 1;
  4981.   do
  4982.     {
  4983.       if (exp & nexp)
  4984.     emul (p, xt, xt);
  4985.       p -= NE;
  4986.       exp = exp + exp;
  4987.     }
  4988.   while (exp <= MAXP);
  4989.  
  4990.   emovi (xt, tt);
  4991.   if (esign < 0)
  4992.     {
  4993.       lexp -= tt[E];
  4994.       k = edivm (tt, yy);
  4995.       lexp += EXONE;
  4996.     }
  4997.   else
  4998.     {
  4999.       lexp += tt[E];
  5000.       k = emulm (tt, yy);
  5001.       lexp -= EXONE - 1;
  5002.     }
  5003.  
  5004.  expdon:
  5005.  
  5006.   /* Round and convert directly to the destination type */
  5007.   if (oprec == 53)
  5008.     lexp -= EXONE - 0x3ff;
  5009. #ifdef IBM
  5010.   else if (oprec == 24 || oprec == 56)
  5011.     lexp -= EXONE - (0x41 << 2);
  5012. #else
  5013.   else if (oprec == 24)
  5014.     lexp -= EXONE - 0177;
  5015. #endif
  5016. #ifdef DEC
  5017.   else if (oprec == 56)
  5018.     lexp -= EXONE - 0201;
  5019. #endif
  5020.   rndprc = oprec;
  5021.   emdnorm (yy, k, 0, lexp, 64);
  5022.  
  5023.  aexit:
  5024.  
  5025.   rndprc = rndsav;
  5026.   yy[0] = nsign;
  5027.   switch (oprec)
  5028.     {
  5029. #ifdef DEC
  5030.     case 56:
  5031.       todec (yy, y);        /* see etodec.c */
  5032.       break;
  5033. #endif
  5034. #ifdef IBM
  5035.     case 56:
  5036.       toibm (yy, y, DFmode);
  5037.       break;
  5038. #endif
  5039.     case 53:
  5040.       toe53 (yy, y);
  5041.       break;
  5042.     case 24:
  5043.       toe24 (yy, y);
  5044.       break;
  5045.     case 64:
  5046.       toe64 (yy, y);
  5047.       break;
  5048.     case 113:
  5049.       toe113 (yy, y);
  5050.       break;
  5051.     case NBITS:
  5052.       emovo (yy, y);
  5053.       break;
  5054.     }
  5055. }
  5056.  
  5057.  
  5058.  
  5059. /* Return Y = largest integer not greater than X (truncated toward minus
  5060.    infinity).  */
  5061.  
  5062. static unsigned EMUSHORT bmask[] =
  5063. {
  5064.   0xffff,
  5065.   0xfffe,
  5066.   0xfffc,
  5067.   0xfff8,
  5068.   0xfff0,
  5069.   0xffe0,
  5070.   0xffc0,
  5071.   0xff80,
  5072.   0xff00,
  5073.   0xfe00,
  5074.   0xfc00,
  5075.   0xf800,
  5076.   0xf000,
  5077.   0xe000,
  5078.   0xc000,
  5079.   0x8000,
  5080.   0x0000,
  5081. };
  5082.  
  5083. static void 
  5084. efloor (x, y)
  5085.      unsigned EMUSHORT x[], y[];
  5086. {
  5087.   register unsigned EMUSHORT *p;
  5088.   int e, expon, i;
  5089.   unsigned EMUSHORT f[NE];
  5090.  
  5091.   emov (x, f);            /* leave in external format */
  5092.   expon = (int) f[NE - 1];
  5093.   e = (expon & 0x7fff) - (EXONE - 1);
  5094.   if (e <= 0)
  5095.     {
  5096.       eclear (y);
  5097.       goto isitneg;
  5098.     }
  5099.   /* number of bits to clear out */
  5100.   e = NBITS - e;
  5101.   emov (f, y);
  5102.   if (e <= 0)
  5103.     return;
  5104.  
  5105.   p = &y[0];
  5106.   while (e >= 16)
  5107.     {
  5108.       *p++ = 0;
  5109.       e -= 16;
  5110.     }
  5111.   /* clear the remaining bits */
  5112.   *p &= bmask[e];
  5113.   /* truncate negatives toward minus infinity */
  5114.  isitneg:
  5115.  
  5116.   if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
  5117.     {
  5118.       for (i = 0; i < NE - 1; i++)
  5119.     {
  5120.       if (f[i] != y[i])
  5121.         {
  5122.           esub (eone, y, y);
  5123.           break;
  5124.         }
  5125.     }
  5126.     }
  5127. }
  5128.  
  5129.  
  5130. /* Return S and EXP such that  S * 2^EXP = X and .5 <= S < 1.
  5131.    For example, 1.1 = 0.55 * 2^1.  */
  5132.  
  5133. static void 
  5134. efrexp (x, exp, s)
  5135.      unsigned EMUSHORT x[];
  5136.      int *exp;
  5137.      unsigned EMUSHORT s[];
  5138. {
  5139.   unsigned EMUSHORT xi[NI];
  5140.   EMULONG li;
  5141.  
  5142.   emovi (x, xi);
  5143.   /*  Handle denormalized numbers properly using long integer exponent.  */
  5144.   li = (EMULONG) ((EMUSHORT) xi[1]);
  5145.  
  5146.   if (li == 0)
  5147.     {
  5148.       li -= enormlz (xi);
  5149.     }
  5150.   xi[1] = 0x3ffe;
  5151.   emovo (xi, s);
  5152.   *exp = (int) (li - 0x3ffe);
  5153. }
  5154.  
  5155. /* Return e type Y = X * 2^PWR2.  */
  5156.  
  5157. static void 
  5158. eldexp (x, pwr2, y)
  5159.      unsigned EMUSHORT x[];
  5160.      int pwr2;
  5161.      unsigned EMUSHORT y[];
  5162. {
  5163.   unsigned EMUSHORT xi[NI];
  5164.   EMULONG li;
  5165.   int i;
  5166.  
  5167.   emovi (x, xi);
  5168.   li = xi[1];
  5169.   li += pwr2;
  5170.   i = 0;
  5171.   emdnorm (xi, i, i, li, 64);
  5172.   emovo (xi, y);
  5173. }
  5174.  
  5175.  
  5176. /* C = remainder after dividing B by A, all e type values.
  5177.    Least significant integer quotient bits left in EQUOT.  */
  5178.  
  5179. static void 
  5180. eremain (a, b, c)
  5181.      unsigned EMUSHORT a[], b[], c[];
  5182. {
  5183.   unsigned EMUSHORT den[NI], num[NI];
  5184.  
  5185. #ifdef NANS
  5186.   if (eisinf (b)
  5187.       || (ecmp (a, ezero) == 0)
  5188.       || eisnan (a)
  5189.       || eisnan (b))
  5190.     {
  5191.       enan (c, 0);
  5192.       return;
  5193.     }
  5194. #endif
  5195.   if (ecmp (a, ezero) == 0)
  5196.     {
  5197.       mtherr ("eremain", SING);
  5198.       eclear (c);
  5199.       return;
  5200.     }
  5201.   emovi (a, den);
  5202.   emovi (b, num);
  5203.   eiremain (den, num);
  5204.   /* Sign of remainder = sign of quotient */
  5205.   if (a[0] == b[0])
  5206.     num[0] = 0;
  5207.   else
  5208.     num[0] = 0xffff;
  5209.   emovo (num, c);
  5210. }
  5211.  
  5212. /*  Return quotient of exploded e-types NUM / DEN in EQUOT,
  5213.     remainder in NUM.  */
  5214.  
  5215. static void 
  5216. eiremain (den, num)
  5217.      unsigned EMUSHORT den[], num[];
  5218. {
  5219.   EMULONG ld, ln;
  5220.   unsigned EMUSHORT j;
  5221.  
  5222.   ld = den[E];
  5223.   ld -= enormlz (den);
  5224.   ln = num[E];
  5225.   ln -= enormlz (num);
  5226.   ecleaz (equot);
  5227.   while (ln >= ld)
  5228.     {
  5229.       if (ecmpm (den, num) <= 0)
  5230.     {
  5231.       esubm (den, num);
  5232.       j = 1;
  5233.     }
  5234.       else
  5235.       j = 0;
  5236.       eshup1 (equot);
  5237.       equot[NI - 1] |= j;
  5238.       eshup1 (num);
  5239.       ln -= 1;
  5240.     }
  5241.   emdnorm (num, 0, 0, ln, 0);
  5242. }
  5243.  
  5244. /* Report an error condition CODE encountered in function NAME.
  5245.    CODE is one of the following:
  5246.  
  5247.     Mnemonic        Value          Significance
  5248.  
  5249.      DOMAIN            1       argument domain error
  5250.      SING              2       function singularity
  5251.      OVERFLOW          3       overflow range error
  5252.      UNDERFLOW         4       underflow range error
  5253.      TLOSS             5       total loss of precision
  5254.      PLOSS             6       partial loss of precision
  5255.      INVALID           7       NaN - producing operation
  5256.      EDOM             33       Unix domain error code
  5257.      ERANGE           34       Unix range error code
  5258.  
  5259.    The order of appearance of the following messages is bound to the
  5260.    error codes defined above.  */
  5261.  
  5262. #define NMSGS 8
  5263. static char *ermsg[NMSGS] =
  5264. {
  5265.   "unknown",            /* error code 0 */
  5266.   "domain",            /* error code 1 */
  5267.   "singularity",        /* et seq.      */
  5268.   "overflow",
  5269.   "underflow",
  5270.   "total loss of precision",
  5271.   "partial loss of precision",
  5272.   "invalid operation"
  5273. };
  5274.  
  5275. int merror = 0;
  5276. extern int merror;
  5277.  
  5278. static void 
  5279. mtherr (name, code)
  5280.      char *name;
  5281.      int code;
  5282. {
  5283.   char errstr[80];
  5284.  
  5285.   /* The string passed by the calling program is supposed to be the
  5286.      name of the function in which the error occurred.
  5287.      The code argument selects which error message string will be printed.  */
  5288.  
  5289.   if ((code <= 0) || (code >= NMSGS))
  5290.     code = 0;
  5291.   sprintf (errstr, " %s %s error", name, ermsg[code]);
  5292.   if (extra_warnings)
  5293.     warning (errstr);
  5294.   /* Set global error message word */
  5295.   merror = code + 1;
  5296. }
  5297.  
  5298. #ifdef DEC
  5299. /* Convert DEC double precision D to e type E.  */
  5300.  
  5301. static void 
  5302. dectoe (d, e)
  5303.      unsigned EMUSHORT *d;
  5304.      unsigned EMUSHORT *e;
  5305. {
  5306.   unsigned EMUSHORT y[NI];
  5307.   register unsigned EMUSHORT r, *p;
  5308.  
  5309.   ecleaz (y);            /* start with a zero */
  5310.   p = y;            /* point to our number */
  5311.   r = *d;            /* get DEC exponent word */
  5312.   if (*d & (unsigned int) 0x8000)
  5313.     *p = 0xffff;        /* fill in our sign */
  5314.   ++p;                /* bump pointer to our exponent word */
  5315.   r &= 0x7fff;            /* strip the sign bit */
  5316.   if (r == 0)            /* answer = 0 if high order DEC word = 0 */
  5317.     goto done;
  5318.  
  5319.  
  5320.   r >>= 7;            /* shift exponent word down 7 bits */
  5321.   r += EXONE - 0201;        /* subtract DEC exponent offset */
  5322.   /* add our e type exponent offset */
  5323.   *p++ = r;            /* to form our exponent */
  5324.  
  5325.   r = *d++;            /* now do the high order mantissa */
  5326.   r &= 0177;            /* strip off the DEC exponent and sign bits */
  5327.   r |= 0200;            /* the DEC understood high order mantissa bit */
  5328.   *p++ = r;            /* put result in our high guard word */
  5329.  
  5330.   *p++ = *d++;            /* fill in the rest of our mantissa */
  5331.   *p++ = *d++;
  5332.   *p = *d;
  5333.  
  5334.   eshdn8 (y);            /* shift our mantissa down 8 bits */
  5335.  done:
  5336.   emovo (y, e);
  5337. }
  5338.  
  5339. /* Convert e type X to DEC double precision D.  */
  5340.  
  5341. static void 
  5342. etodec (x, d)
  5343.      unsigned EMUSHORT *x, *d;
  5344. {
  5345.   unsigned EMUSHORT xi[NI];
  5346.   EMULONG exp;
  5347.   int rndsav;
  5348.  
  5349.   emovi (x, xi);
  5350.   /* Adjust exponent for offsets.  */
  5351.   exp = (EMULONG) xi[E] - (EXONE - 0201);
  5352.   /* Round off to nearest or even.  */
  5353.   rndsav = rndprc;
  5354.   rndprc = 56;
  5355.   emdnorm (xi, 0, 0, exp, 64);
  5356.   rndprc = rndsav;
  5357.   todec (xi, d);
  5358. }
  5359.  
  5360. /* Convert exploded e-type X, that has already been rounded to
  5361.    56-bit precision, to DEC format double Y.  */
  5362.  
  5363. static void 
  5364. todec (x, y)
  5365.      unsigned EMUSHORT *x, *y;
  5366. {
  5367.   unsigned EMUSHORT i;
  5368.   unsigned EMUSHORT *p;
  5369.  
  5370.   p = x;
  5371.   *y = 0;
  5372.   if (*p++)
  5373.     *y = 0100000;
  5374.   i = *p++;
  5375.   if (i == 0)
  5376.     {
  5377.       *y++ = 0;
  5378.       *y++ = 0;
  5379.       *y++ = 0;
  5380.       *y++ = 0;
  5381.       return;
  5382.     }
  5383.   if (i > 0377)
  5384.     {
  5385.       *y++ |= 077777;
  5386.       *y++ = 0xffff;
  5387.       *y++ = 0xffff;
  5388.       *y++ = 0xffff;
  5389. #ifdef ERANGE
  5390.       errno = ERANGE;
  5391. #endif
  5392.       return;
  5393.     }
  5394.   i &= 0377;
  5395.   i <<= 7;
  5396.   eshup8 (x);
  5397.   x[M] &= 0177;
  5398.   i |= x[M];
  5399.   *y++ |= i;
  5400.   *y++ = x[M + 1];
  5401.   *y++ = x[M + 2];
  5402.   *y++ = x[M + 3];
  5403. }
  5404. #endif /* DEC */
  5405.  
  5406. #ifdef IBM
  5407. /* Convert IBM single/double precision to e type.  */
  5408.  
  5409. static void 
  5410. ibmtoe (d, e, mode)
  5411.      unsigned EMUSHORT *d;
  5412.      unsigned EMUSHORT *e;
  5413.      enum machine_mode mode;
  5414. {
  5415.   unsigned EMUSHORT y[NI];
  5416.   register unsigned EMUSHORT r, *p;
  5417.   int rndsav;
  5418.  
  5419.   ecleaz (y);            /* start with a zero */
  5420.   p = y;            /* point to our number */
  5421.   r = *d;            /* get IBM exponent word */
  5422.   if (*d & (unsigned int) 0x8000)
  5423.     *p = 0xffff;        /* fill in our sign */
  5424.   ++p;                /* bump pointer to our exponent word */
  5425.   r &= 0x7f00;            /* strip the sign bit */
  5426.   r >>= 6;            /* shift exponent word down 6 bits */
  5427.                 /* in fact shift by 8 right and 2 left */
  5428.   r += EXONE - (0x41 << 2);    /* subtract IBM exponent offset */
  5429.                   /* add our e type exponent offset */
  5430.   *p++ = r;            /* to form our exponent */
  5431.  
  5432.   *p++ = *d++ & 0xff;        /* now do the high order mantissa */
  5433.                 /* strip off the IBM exponent and sign bits */
  5434.   if (mode != SFmode)        /* there are only 2 words in SFmode */
  5435.     {
  5436.       *p++ = *d++;        /* fill in the rest of our mantissa */
  5437.       *p++ = *d++;
  5438.     }
  5439.   *p = *d;
  5440.  
  5441.   if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
  5442.     y[0] = y[E] = 0;
  5443.   else
  5444.     y[E] -= 5 + enormlz (y);    /* now normalise the mantissa */
  5445.                   /* handle change in RADIX */
  5446.   emovo (y, e);
  5447. }
  5448.  
  5449.  
  5450.  
  5451. /* Convert e type to IBM single/double precision.  */
  5452.  
  5453. static void 
  5454. etoibm (x, d, mode)
  5455.      unsigned EMUSHORT *x, *d;
  5456.      enum machine_mode mode;
  5457. {
  5458.   unsigned EMUSHORT xi[NI];
  5459.   EMULONG exp;
  5460.   int rndsav;
  5461.  
  5462.   emovi (x, xi);
  5463.   exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2));    /* adjust exponent for offsets */
  5464.                             /* round off to nearest or even */
  5465.   rndsav = rndprc;
  5466.   rndprc = 56;
  5467.   emdnorm (xi, 0, 0, exp, 64);
  5468.   rndprc = rndsav;
  5469.   toibm (xi, d, mode);
  5470. }
  5471.  
  5472. static void 
  5473. toibm (x, y, mode)
  5474.      unsigned EMUSHORT *x, *y;
  5475.      enum machine_mode mode;
  5476. {
  5477.   unsigned EMUSHORT i;
  5478.   unsigned EMUSHORT *p;
  5479.   int r;
  5480.  
  5481.   p = x;
  5482.   *y = 0;
  5483.   if (*p++)
  5484.     *y = 0x8000;
  5485.   i = *p++;
  5486.   if (i == 0)
  5487.     {
  5488.       *y++ = 0;
  5489.       *y++ = 0;
  5490.       if (mode != SFmode)
  5491.     {
  5492.       *y++ = 0;
  5493.       *y++ = 0;
  5494.     }
  5495.       return;
  5496.     }
  5497.   r = i & 0x3;
  5498.   i >>= 2;
  5499.   if (i > 0x7f)
  5500.     {
  5501.       *y++ |= 0x7fff;
  5502.       *y++ = 0xffff;
  5503.       if (mode != SFmode)
  5504.     {
  5505.       *y++ = 0xffff;
  5506.       *y++ = 0xffff;
  5507.     }
  5508. #ifdef ERANGE
  5509.       errno = ERANGE;
  5510. #endif
  5511.       return;
  5512.     }
  5513.   i &= 0x7f;
  5514.   *y |= (i << 8);
  5515.   eshift (x, r + 5);
  5516.   *y++ |= x[M];
  5517.   *y++ = x[M + 1];
  5518.   if (mode != SFmode)
  5519.     {
  5520.       *y++ = x[M + 2];
  5521.       *y++ = x[M + 3];
  5522.     }
  5523. }
  5524. #endif /* IBM */
  5525.  
  5526. /* Output a binary NaN bit pattern in the target machine's format.  */
  5527.  
  5528. /* If special NaN bit patterns are required, define them in tm.h
  5529.    as arrays of unsigned 16-bit shorts.  Otherwise, use the default
  5530.    patterns here. */
  5531. #ifdef TFMODE_NAN
  5532. TFMODE_NAN;
  5533. #else
  5534. #ifdef IEEE
  5535. unsigned EMUSHORT TFbignan[8] =
  5536.  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
  5537. unsigned EMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
  5538. #endif
  5539. #endif
  5540.  
  5541. #ifdef XFMODE_NAN
  5542. XFMODE_NAN;
  5543. #else
  5544. #ifdef IEEE
  5545. unsigned EMUSHORT XFbignan[6] =
  5546.  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
  5547. unsigned EMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
  5548. #endif
  5549. #endif
  5550.  
  5551. #ifdef DFMODE_NAN
  5552. DFMODE_NAN;
  5553. #else
  5554. #ifdef IEEE
  5555. unsigned EMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
  5556. unsigned EMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
  5557. #endif
  5558. #endif
  5559.  
  5560. #ifdef SFMODE_NAN
  5561. SFMODE_NAN;
  5562. #else
  5563. #ifdef IEEE
  5564. unsigned EMUSHORT SFbignan[2] = {0x7fff, 0xffff};
  5565. unsigned EMUSHORT SFlittlenan[2] = {0, 0xffc0};
  5566. #endif
  5567. #endif
  5568.  
  5569.  
  5570. static void
  5571. make_nan (nan, sign, mode)
  5572.      unsigned EMUSHORT *nan;
  5573.      int sign;
  5574.      enum machine_mode mode;
  5575. {
  5576.   int n;
  5577.   unsigned EMUSHORT *p;
  5578.  
  5579.   switch (mode)
  5580.     {
  5581. /* Possibly the `reserved operand' patterns on a VAX can be
  5582.    used like NaN's, but probably not in the same way as IEEE. */
  5583. #if !defined(DEC) && !defined(IBM)
  5584.     case TFmode:
  5585.       n = 8;
  5586.       if (REAL_WORDS_BIG_ENDIAN)
  5587.     p = TFbignan;
  5588.       else
  5589.     p = TFlittlenan;
  5590.       break;
  5591.     case XFmode:
  5592.       n = 6;
  5593.       if (REAL_WORDS_BIG_ENDIAN)
  5594.     p = XFbignan;
  5595.       else
  5596.     p = XFlittlenan;
  5597.       break;
  5598.     case DFmode:
  5599.       n = 4;
  5600.       if (REAL_WORDS_BIG_ENDIAN)
  5601.     p = DFbignan;
  5602.       else
  5603.     p = DFlittlenan;
  5604.       break;
  5605.     case HFmode:
  5606.     case SFmode:
  5607.       n = 2;
  5608.       if (REAL_WORDS_BIG_ENDIAN)
  5609.     p = SFbignan;
  5610.       else
  5611.     p = SFlittlenan;
  5612.       break;
  5613. #endif
  5614.     default:
  5615.       abort ();
  5616.     }
  5617.   if (REAL_WORDS_BIG_ENDIAN)
  5618.     *nan++ = (sign << 15) | *p++;
  5619.   while (--n != 0)
  5620.     *nan++ = *p++;
  5621.   if (! REAL_WORDS_BIG_ENDIAN)
  5622.     *nan = (sign << 15) | *p;
  5623. }
  5624.  
  5625. /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
  5626.    This is the inverse of the function `etarsingle' invoked by
  5627.    REAL_VALUE_TO_TARGET_SINGLE.  */
  5628.  
  5629. REAL_VALUE_TYPE
  5630. ereal_from_float (f)
  5631.      HOST_WIDE_INT f;
  5632. {
  5633.   REAL_VALUE_TYPE r;
  5634.   unsigned EMUSHORT s[2];
  5635.   unsigned EMUSHORT e[NE];
  5636.  
  5637.   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
  5638.    This is the inverse operation to what the function `endian' does.  */
  5639.   if (REAL_WORDS_BIG_ENDIAN)
  5640.     {
  5641.       s[0] = (unsigned EMUSHORT) (f >> 16);
  5642.       s[1] = (unsigned EMUSHORT) f;
  5643.     }
  5644.   else
  5645.     {
  5646.       s[0] = (unsigned EMUSHORT) f;
  5647.       s[1] = (unsigned EMUSHORT) (f >> 16);
  5648.     }
  5649.   /* Convert and promote the target float to E-type. */
  5650.   e24toe (s, e);
  5651.   /* Output E-type to REAL_VALUE_TYPE. */
  5652.   PUT_REAL (e, &r);
  5653.   return r;
  5654. }
  5655.  
  5656.  
  5657. /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
  5658.    This is the inverse of the function `etardouble' invoked by
  5659.    REAL_VALUE_TO_TARGET_DOUBLE.
  5660.  
  5661.    The DFmode is stored as an array of HOST_WIDE_INT in the target's
  5662.    data format, with no holes in the bit packing.  The first element
  5663.    of the input array holds the bits that would come first in the
  5664.    target computer's memory.  */
  5665.  
  5666. REAL_VALUE_TYPE
  5667. ereal_from_double (d)
  5668.      HOST_WIDE_INT d[];
  5669. {
  5670.   REAL_VALUE_TYPE r;
  5671.   unsigned EMUSHORT s[4];
  5672.   unsigned EMUSHORT e[NE];
  5673.  
  5674.   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
  5675.   if (REAL_WORDS_BIG_ENDIAN)
  5676.     {
  5677.       s[0] = (unsigned EMUSHORT) (d[0] >> 16);
  5678.       s[1] = (unsigned EMUSHORT) d[0];
  5679. #if HOST_BITS_PER_WIDE_INT == 32
  5680.       s[2] = (unsigned EMUSHORT) (d[1] >> 16);
  5681.       s[3] = (unsigned EMUSHORT) d[1];
  5682. #else
  5683.       /* In this case the entire target double is contained in the
  5684.      first array element.  The second element of the input is
  5685.      ignored.  */
  5686.       s[2] = (unsigned EMUSHORT) (d[0] >> 48);
  5687.       s[3] = (unsigned EMUSHORT) (d[0] >> 32);
  5688. #endif
  5689.     }
  5690.   else
  5691.     {
  5692.       /* Target float words are little-endian.  */
  5693.       s[0] = (unsigned EMUSHORT) d[0];
  5694.       s[1] = (unsigned EMUSHORT) (d[0] >> 16);
  5695. #if HOST_BITS_PER_WIDE_INT == 32
  5696.       s[2] = (unsigned EMUSHORT) d[1];
  5697.       s[3] = (unsigned EMUSHORT) (d[1] >> 16);
  5698. #else
  5699.       s[2] = (unsigned EMUSHORT) (d[0] >> 32);
  5700.       s[3] = (unsigned EMUSHORT) (d[0] >> 48);
  5701. #endif
  5702.     }
  5703.   /* Convert target double to E-type. */
  5704.   e53toe (s, e);
  5705.   /* Output E-type to REAL_VALUE_TYPE. */
  5706.   PUT_REAL (e, &r);
  5707.   return r;
  5708. }
  5709.  
  5710.  
  5711. /* Convert target computer unsigned 64-bit integer to e-type.
  5712.    The endian-ness of DImode follows the convention for integers,
  5713.    so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN.  */
  5714.  
  5715. static void
  5716. uditoe (di, e)
  5717.      unsigned EMUSHORT *di;  /* Address of the 64-bit int. */
  5718.      unsigned EMUSHORT *e;
  5719. {
  5720.   unsigned EMUSHORT yi[NI];
  5721.   int k;
  5722.  
  5723.   ecleaz (yi);
  5724.   if (WORDS_BIG_ENDIAN)
  5725.     {
  5726.       for (k = M; k < M + 4; k++)
  5727.     yi[k] = *di++;
  5728.     }
  5729.   else
  5730.     {
  5731.       for (k = M + 3; k >= M; k--)
  5732.     yi[k] = *di++;
  5733.     }
  5734.   yi[E] = EXONE + 47;    /* exponent if normalize shift count were 0 */
  5735.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  5736.     ecleaz (yi);        /* it was zero */
  5737.   else
  5738.     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
  5739.   emovo (yi, e);
  5740. }
  5741.  
  5742. /* Convert target computer signed 64-bit integer to e-type. */
  5743.  
  5744. static void
  5745. ditoe (di, e)
  5746.      unsigned EMUSHORT *di;  /* Address of the 64-bit int. */
  5747.      unsigned EMUSHORT *e;
  5748. {
  5749.   unsigned EMULONG acc;
  5750.   unsigned EMUSHORT yi[NI];
  5751.   unsigned EMUSHORT carry;
  5752.   int k, sign;
  5753.  
  5754.   ecleaz (yi);
  5755.   if (WORDS_BIG_ENDIAN)
  5756.     {
  5757.       for (k = M; k < M + 4; k++)
  5758.     yi[k] = *di++;
  5759.     }
  5760.   else
  5761.     {
  5762.       for (k = M + 3; k >= M; k--)
  5763.     yi[k] = *di++;
  5764.     }
  5765.   /* Take absolute value */
  5766.   sign = 0;
  5767.   if (yi[M] & 0x8000)
  5768.     {
  5769.       sign = 1;
  5770.       carry = 0;
  5771.       for (k = M + 3; k >= M; k--)
  5772.     {
  5773.       acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
  5774.       yi[k] = acc;
  5775.       carry = 0;
  5776.       if (acc & 0x10000)
  5777.         carry = 1;
  5778.     }
  5779.     }
  5780.   yi[E] = EXONE + 47;    /* exponent if normalize shift count were 0 */
  5781.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  5782.     ecleaz (yi);        /* it was zero */
  5783.   else
  5784.     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
  5785.   emovo (yi, e);
  5786.   if (sign)
  5787.     eneg (e);
  5788. }
  5789.  
  5790.  
  5791. /* Convert e-type to unsigned 64-bit int. */
  5792.  
  5793. static void 
  5794. etoudi (x, i)
  5795.      unsigned EMUSHORT *x;
  5796.      unsigned EMUSHORT *i;
  5797. {
  5798.   unsigned EMUSHORT xi[NI];
  5799.   int j, k;
  5800.  
  5801.   emovi (x, xi);
  5802.   if (xi[0])
  5803.     {
  5804.       xi[M] = 0;
  5805.       goto noshift;
  5806.     }
  5807.   k = (int) xi[E] - (EXONE - 1);
  5808.   if (k <= 0)
  5809.     {
  5810.       for (j = 0; j < 4; j++)
  5811.     *i++ = 0;
  5812.       return;
  5813.     }
  5814.   if (k > 64)
  5815.     {
  5816.       for (j = 0; j < 4; j++)
  5817.     *i++ = 0xffff;
  5818.       if (extra_warnings)
  5819.     warning ("overflow on truncation to integer");
  5820.       return;
  5821.     }
  5822.   if (k > 16)
  5823.     {
  5824.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  5825.      then shift up by 16's.  */
  5826.       j = k - ((k >> 4) << 4);
  5827.       if (j == 0)
  5828.     j = 16;
  5829.       eshift (xi, j);
  5830.       if (WORDS_BIG_ENDIAN)
  5831.     *i++ = xi[M];
  5832.       else
  5833.     {
  5834.       i += 3;
  5835.       *i-- = xi[M];
  5836.     }
  5837.       k -= j;
  5838.       do
  5839.     {
  5840.       eshup6 (xi);
  5841.       if (WORDS_BIG_ENDIAN)
  5842.         *i++ = xi[M];
  5843.       else
  5844.         *i-- = xi[M];
  5845.     }
  5846.       while ((k -= 16) > 0);
  5847.     }
  5848.   else
  5849.     {
  5850.         /* shift not more than 16 bits */
  5851.       eshift (xi, k);
  5852.  
  5853. noshift:
  5854.  
  5855.       if (WORDS_BIG_ENDIAN)
  5856.     {
  5857.       i += 3;
  5858.       *i-- = xi[M];
  5859.       *i-- = 0;
  5860.       *i-- = 0;
  5861.       *i = 0;
  5862.     }
  5863.       else
  5864.     {
  5865.       *i++ = xi[M];
  5866.       *i++ = 0;
  5867.       *i++ = 0;
  5868.       *i = 0;
  5869.     }
  5870.     }
  5871. }
  5872.  
  5873.  
  5874. /* Convert e-type to signed 64-bit int. */
  5875.  
  5876. static void 
  5877. etodi (x, i)
  5878.      unsigned EMUSHORT *x;
  5879.      unsigned EMUSHORT *i;
  5880. {
  5881.   unsigned EMULONG acc;
  5882.   unsigned EMUSHORT xi[NI];
  5883.   unsigned EMUSHORT carry;
  5884.   unsigned EMUSHORT *isave;
  5885.   int j, k;
  5886.  
  5887.   emovi (x, xi);
  5888.   k = (int) xi[E] - (EXONE - 1);
  5889.   if (k <= 0)
  5890.     {
  5891.       for (j = 0; j < 4; j++)
  5892.     *i++ = 0;
  5893.       return;
  5894.     }
  5895.   if (k > 64)
  5896.     {
  5897.       for (j = 0; j < 4; j++)
  5898.     *i++ = 0xffff;
  5899.       if (extra_warnings)
  5900.     warning ("overflow on truncation to integer");
  5901.       return;
  5902.     }
  5903.   isave = i;
  5904.   if (k > 16)
  5905.     {
  5906.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  5907.      then shift up by 16's.  */
  5908.       j = k - ((k >> 4) << 4);
  5909.       if (j == 0)
  5910.     j = 16;
  5911.       eshift (xi, j);
  5912.       if (WORDS_BIG_ENDIAN)
  5913.     *i++ = xi[M];
  5914.       else
  5915.     {
  5916.       i += 3;
  5917.       *i-- = xi[M];
  5918.     }
  5919.       k -= j;
  5920.       do
  5921.     {
  5922.       eshup6 (xi);
  5923.       if (WORDS_BIG_ENDIAN)
  5924.         *i++ = xi[M];
  5925.       else
  5926.         *i-- = xi[M];
  5927.     }
  5928.       while ((k -= 16) > 0);
  5929.     }
  5930.   else
  5931.     {
  5932.         /* shift not more than 16 bits */
  5933.       eshift (xi, k);
  5934.  
  5935.       if (WORDS_BIG_ENDIAN)
  5936.     {
  5937.       i += 3;
  5938.       *i = xi[M];
  5939.       *i-- = 0;
  5940.       *i-- = 0;
  5941.       *i = 0;
  5942.     }
  5943.       else
  5944.     {
  5945.       *i++ = xi[M];
  5946.       *i++ = 0;
  5947.       *i++ = 0;
  5948.       *i = 0;
  5949.     }
  5950.     }
  5951.   /* Negate if negative */
  5952.   if (xi[0])
  5953.     {
  5954.       carry = 0;
  5955.       if (WORDS_BIG_ENDIAN)
  5956.     isave += 3;
  5957.       for (k = 0; k < 4; k++)
  5958.     {
  5959.       acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
  5960.       if (WORDS_BIG_ENDIAN)
  5961.         *isave-- = acc;
  5962.       else
  5963.         *isave++ = acc;
  5964.       carry = 0;
  5965.       if (acc & 0x10000)
  5966.         carry = 1;
  5967.     }
  5968.     }
  5969. }
  5970.  
  5971.  
  5972. /* Longhand square root routine. */
  5973.  
  5974.  
  5975. static int esqinited = 0;
  5976. static unsigned short sqrndbit[NI];
  5977.  
  5978. static void 
  5979. esqrt (x, y)
  5980.      unsigned EMUSHORT *x, *y;
  5981. {
  5982.   unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
  5983.   EMULONG m, exp;
  5984.   int i, j, k, n, nlups;
  5985.  
  5986.   if (esqinited == 0)
  5987.     {
  5988.       ecleaz (sqrndbit);
  5989.       sqrndbit[NI - 2] = 1;
  5990.       esqinited = 1;
  5991.     }
  5992.   /* Check for arg <= 0 */
  5993.   i = ecmp (x, ezero);
  5994.   if (i <= 0)
  5995.     {
  5996.       if (i == -1)
  5997.     {
  5998.       mtherr ("esqrt", DOMAIN);
  5999.       eclear (y);
  6000.     }
  6001.       else
  6002.     emov (x, y);
  6003.       return;
  6004.     }
  6005.  
  6006. #ifdef INFINITY
  6007.   if (eisinf (x))
  6008.     {
  6009.       eclear (y);
  6010.       einfin (y);
  6011.       return;
  6012.     }
  6013. #endif
  6014.   /* Bring in the arg and renormalize if it is denormal. */
  6015.   emovi (x, xx);
  6016.   m = (EMULONG) xx[1];        /* local long word exponent */
  6017.   if (m == 0)
  6018.     m -= enormlz (xx);
  6019.  
  6020.   /* Divide exponent by 2 */
  6021.   m -= 0x3ffe;
  6022.   exp = (unsigned short) ((m / 2) + 0x3ffe);
  6023.  
  6024.   /* Adjust if exponent odd */
  6025.   if ((m & 1) != 0)
  6026.     {
  6027.       if (m > 0)
  6028.     exp += 1;
  6029.       eshdn1 (xx);
  6030.     }
  6031.  
  6032.   ecleaz (sq);
  6033.   ecleaz (num);
  6034.   n = 8;            /* get 8 bits of result per inner loop */
  6035.   nlups = rndprc;
  6036.   j = 0;
  6037.  
  6038.   while (nlups > 0)
  6039.     {
  6040.       /* bring in next word of arg */
  6041.       if (j < NE)
  6042.     num[NI - 1] = xx[j + 3];
  6043.       /* Do additional bit on last outer loop, for roundoff. */
  6044.       if (nlups <= 8)
  6045.     n = nlups + 1;
  6046.       for (i = 0; i < n; i++)
  6047.     {
  6048.       /* Next 2 bits of arg */
  6049.       eshup1 (num);
  6050.       eshup1 (num);
  6051.       /* Shift up answer */
  6052.       eshup1 (sq);
  6053.       /* Make trial divisor */
  6054.       for (k = 0; k < NI; k++)
  6055.         temp[k] = sq[k];
  6056.       eshup1 (temp);
  6057.       eaddm (sqrndbit, temp);
  6058.       /* Subtract and insert answer bit if it goes in */
  6059.       if (ecmpm (temp, num) <= 0)
  6060.         {
  6061.           esubm (temp, num);
  6062.           sq[NI - 2] |= 1;
  6063.         }
  6064.     }
  6065.       nlups -= n;
  6066.       j += 1;
  6067.     }
  6068.  
  6069.   /* Adjust for extra, roundoff loop done. */
  6070.   exp += (NBITS - 1) - rndprc;
  6071.  
  6072.   /* Sticky bit = 1 if the remainder is nonzero. */
  6073.   k = 0;
  6074.   for (i = 3; i < NI; i++)
  6075.     k |= (int) num[i];
  6076.  
  6077.   /* Renormalize and round off. */
  6078.   emdnorm (sq, k, 0, exp, 64);
  6079.   emovo (sq, y);
  6080. }
  6081. #endif /* EMU_NON_COMPILE not defined */
  6082.  
  6083. /* Return the binary precision of the significand for a given
  6084.    floating point mode.  The mode can hold an integer value
  6085.    that many bits wide, without losing any bits.  */
  6086.  
  6087. int
  6088. significand_size (mode)
  6089.      enum machine_mode mode;
  6090. {
  6091.  
  6092. switch (mode)
  6093.   {
  6094.   case SFmode:
  6095.     return 24;
  6096.  
  6097.   case DFmode:
  6098. #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
  6099.     return 53;
  6100. #else
  6101. #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
  6102.     return 56;
  6103. #else
  6104. #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  6105.     return 56;
  6106. #else
  6107.     abort ();
  6108. #endif
  6109. #endif
  6110. #endif
  6111.  
  6112.   case XFmode:
  6113.     return 64;
  6114.   case TFmode:
  6115.     return 113;
  6116.  
  6117.   default:
  6118.     abort ();
  6119.   }
  6120. }
  6121.